ON  THE  AUTOMATED  OPTIMAL  DESIGN 
OF 
CONSTRAINED  STRUCTURES 


BY 
JERRY  C.  KORNBUCKLE 


A  DISSERTATION  PRESENTED  TO  THE  GRADUATE  COUNCIL  OF 
THE  LT^JIVERSITY  OF  FLORIDA  IN  PARTIAL 
FULFLLLKENT  OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 


UNIVERSITY  OF  FLORIDA 
1974 


Copyright  by 
Jerry  C.  Hornbuckle 
1974 


DEDICATION 

To  my  grandmother,  Mrs.  Florence  Hornbuckle, 
and  my  wife,  Carolyn.   Without  the  love,  confidence,  and 
patient  understanding  of  Granny  Hornbuckle  and  Carolyn 
my  graduate  studies  would  never  have  been  attempted. 


ACKNOWLEDGMENTS 

To  Dr.  Robert  L.  Sierakowski  and  Dr.  William  H.  Boykin,  Jr., 
for  guiding  my  research  and  for  being  more  than  just  advisors. 

To  Dr.  Gene  W.  Hemp,  Dr.  Ibrahim  K.  Ebcioglu,  and  Dr.  John 
M.  Vance,  for  their  assistance,  support,  and  for  serving  on  my  advisory 
cominittee. 

To  Dr.  Lawrence  E.  Malvern  and  Dr.  Martin  A.  Eisenberg,  for 
always  finding  the  time  to  offer  advice  and  explanations  on  questions 
related  to  solid  mechanics  and  academics. 

To  the  departmental  office  staff  for  their  kind  assistance 
with  administrative  problems  and  clerical  support. 

To  Randell  A.  Crowe,  Charles  D.  Myers,  and  J.  Eric  Schonblom 
for  attentive  discussions  of  many  little  problems  and  for  assistance 
in  preparing  for  the  qualifying  examination. 


IV 


TABLE  OF  CONTENTS 

Page 

ACKNOWLEDGMENTS  iv 

ABSTRACT viii 

CHAPTER 

I    INTRODUCTION   1 

1.0  Survey  Papers   1 

1.1  Historical  Development:   Optimal  Columns  6 

1.2  Historical  Development:   Optimal  Static  Beams   ...  9 

1.3  Historical  Development:   Optimal  Dynamical  Beams  .  .  12 

1.4  Scope  of  the  Dissertation 14 

II    GENERAL  PROBLEIIS  AND  METHODS  IN  STRUCTUR/.L  OPTIMIZATION  .  16 

2.0  Introduction 15 

2.1  Problem  Classification  Criteria   16 

2.1.0  Problem  Classification  Guidelines  18 

2.1.1  Governing  Equations  of  the  System IS 

2.1.2  Constraints 19 

2.1.3  Cost  Functionals 21 

2.2  Methods:   Continuous  Systems  26 

2.2.0  Special  Variational  Methods  27 

2.2.1  Energy  Methods 27 

2.2.2  Pontryagin's  Maximum  Principle   29 

2.2.3  Method  of  Steepest  Ascent/Descent  31 

2.2.4  Transition  Matrix:   Aeroelasticity  Problems  .  31 

2.2.5  Oth.er  Miscellaneous  Methods 32 

2.3  Methods:   Discrete  Systems  33 

2.3.0  Mathematical  Programming   34 

2.3.1  Discrete  Solution  Approximations   35 

2.3.2  Segmentwise-Constant  Approximations  37 

2.3.3  Complex  Structures  with  Frequency 

Constraints 38 

2.3.4  Finite  Element  Approximations  40 

2.4  Closure 41 


TABLE  OF  CONTENTS  (Continued) 

CHAPTER  Page 

III    THEORETICAL  DEVELOPMENT  43 

3.0  Introduction 43 

3.1  Problem  Statement  and  Necessary  Conditions  44 

3.2  Mathematical  Programming:   Gradient 

Projection  Method  .....  48 

3.3  Gradient  Projection  Methods  Applied  to 

the  Maximum  Principle  53 

3.4  Maximum  Principle  Algorithm   60 

3.5  Solution  Methods  63 

IV    CONSTRAINED  DESIGN  OF  A  CANTILEVER  BEAiM 

BENDING  DUE  TO  ITS  0\-m   kTlIGHT 66 

4.0  Introduction 66 

4.1  Problem  Statement   66 

4.2  Structural  System   67 

4.3  Unmodified  Application  of  the  Maximum  Principle   .  .  70 

4.4  Results:   Geometric  Control  Constraints   79 

4.5  Inequality  Stress  Constraints   89 

4.6  Results:   Stress  Constraints  Included   93 

V    CONSTRrMNED  DESIGN  FOR  AN  OPTIMAL  EIGENVALUE  PROBLEM   .  .  101 

5.0  Introduction 101 

5.1  Problem  Statement   101 

5.2  Structural  System   . 102 

5.3  Analysis  of  the  Problem 109 

5.4  Application  of  the  Maxim.um  Principle  .  , 121 

5.5  Results:   Geometric  Control  Constraints   132 

5.6  Inequality  Stress  Constraints   148 

VI    FIKITE  ELEMENT  METHODS  IN  STRUCTUPJ^L  OPTIMIZATION: 

AN  EXAJ-IPLE 155 

6.0  Introduction 155 

6.1  Finite  Element  Problem  Statement  155 

6.2  Mathematical  Programming:   Gradient 

Projection  Method    157 

6.3  Results 162 

VII    COMMENTS  ON  NTOIERICAL  INSTABILITY  IN  THE 

QUASILINEARIZATION  ALGORITHM  173 

7.0  Introduction 173 

7.1  Computer  Program  Convergence  Features   173 

7.2  Num.erical  Instatilibies  for  Cantilever 

Beam  Example 175 

7.3  Numerical  Instabilities  for  Column 

Buckling  Example   183 

vi 


TABLE  OF  CONTENTS  (Continued) 

CHAPTER  Page 

VIII    CONCLUSIONS  AND  RECOMMENDATIONS  186 

8.0  Summary  and  Conclusions 186 

8.1  Recommendations 186 

APPENDIXES 

A    HISTORICAL  DEVELOPMENTS  191 

B    A  SIMPLE  PROOF  OF  THE  KUHN-TUCKER  THEOREM  206 

C    COMPUTER  SUBROUTINE  LISTINGS   214 

BIBLIOGRAPHY   244 

BIOGRAPHICAL  SKETCH  25  4 


VI 1 


Abstract  of  Dissertation  Presented  to  the  Graduate  Council 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the  Requirements 
for  the  Degree  of  Doctor  of  Philosophy 


ON  THE  AUTOMATED  OPTIMAL  DESIGN 
OF  CONSTRAINED  STRUCTURES 

By 

Jerry  C.  Hornbuckle 

August,  1974 


Chairican:   Dr.  William  H.  Boykin,  Jr. 
Co-chairman:   Dr.  Robert  L.  Sierakowski 
Major  Department:   Engineering  Sciences 


Pontryagin's  Maximum  Principle  is  applied  to  the  optimal  design 
of  elastic  structures,  subject  to  both  hard  inequality  constraints  and 
subsidiary  conditions.   By  analyzing  the  maximum  principle  as  a  non- 
linear programming  problem,  an  explicit  formulation  is  derived  for  the 
Lagrangian  multiplier  functions  that  adjoin  the  constraints  to  the  cost 
functional.   With  this  result  the  usual  necessary  conditions  for  opti- 
mality  can  themselves  be  used  directly  in  an  algorithm  for  obtaining 
a  solution. 

A  survey  of  general  methods  and  problems  in  the  optimal  design 
of  elastic  structures  shows  that  there  are  two  general  problem  types 
depending  upon  whether  or  not  the  cost  functional  is  an  eigenvalue. 
An  example  problem  of  each  type  is  included  with  the  solutions  obtained 
by  the  method  of  quasilinearization.   In  the  first  example,  a  minimum 
deflection  beam  problem,  classical  Maximum  Principle  techniques  are 
used.   The  eigenvalue  problem  is  exemplified  by  the  maximization  of  the 
buckling  load  of  a  column  and  uses  the  explicit  multiplier  function 

viii 


formulation  mentioned  above.   Since  the  problem  considered  is  conserva- 
tive, it  is  therefore  described  mathematically  by  a  self-adjoint  system; 
under  this  condition  it  is  shown  that  the  minimum  weight  problem  is 
identical  to  the  maximum  buckling  load  problem. 

In  order  to  demonstrate  the  theory  for  the  programming  techniques 
used,  the  beam  problem  is  also  solved  by  using  a  finite  element  repre- 
sentation of  the  structure.   From  a  comparison  to  the  maximum  principle 
solution  it  is  shown  that  the  form  of  the  optimal  solution  obtained  is 
dependent  upon  the  magnitude  of  the  tolerance  used  with  the  numerical 
solution  scheme.   Furthermore,  it  is  shown  that  convergence  by  the 
quasilinearization  algorithm  is  related  to  the  respective  curvatures  of 
the  initial  guess  and  the  solution. 

Recommendations  for  additional  investigations  pertinent  to  this 
study  are  also  included. 


XX 


CHAPTER  I 


INTRODUCTION 


1.0   Survey  Papers 

It  is  exceedingly  difficult  to  write  a  general  introduction  to 
the  field  of  structural  optimization  for  two  basic  reasons:   (i)  there 
is  no  conventionally  accepted  nomenclature,  and  (ii)  there  is  also  no 
conventionally  accepted  classification  of  problem  types  or  character- 
istics.  In  marked  contrast,  when  one  considers  the  calculus  of  varia- 
tions, "cost  functional,  system  equations,  kinematic  and  natural  bound- 
ary conditions,  adjoint  variables,  Hamiltonian,  etc.,"  all  have  well- 
defined,  universally  accepted  meanings.   Additionally,  there  is  no  con- 
fusion when  speaking  of  the  problem  types  of  Mayer,  Lagrange,  and  Bolza. 
This  conunon  language  and  categorization  of  problems  does  not  exist  in 
structural  optimization. 

Instead,  the  field  tends  to  branch  and  fragment  into  very  special- 
ized sub-disciplines  that  are  oriented  towards  applications.   While  these 
branches  are  related  to  the  general  field,  the  techniques  and  methods  of 
one  branch  can  seldom  be  applied  to  another.   Moreover,  as  a  result  of 
the  tendency  to  an  applications  orientation,  solutions  are  generally  ad 
lioc  and  not  useful  for  other  problems  even  within  the  same  branch.   The 
lack  of  any  definitive  unification  of  the  subject  cannot  be  blamed  either 
on  being  recently  developed  or  in  receiving  too  little  attention. 


This  is  readily  seen  by  considering  the  survey  papers  described  in  the 
follov7ing  paragraphs. 

The  earliest  comprehensive  survey  paper  is  Wasiutynski  and 
Brandt  (1963).   Although  their  excellent  historical  development  is 
dominated  by  Russian  and  Eastern  European  references,  the  authors  do 
include  a  higher  percentage  of  papers  by  Western  authors  than  is 
encountered  in  the  typical  paper  from  Eastern  Europe.   A  more  funda- 
mental criticism  is  that  too  little  is  said  regarding  problem  types  or 
solution  methods. 

Chronologically,  the  next  survey  paper  is  Gerard  (1966)  . 
The  theme  of  this  review  is  aerospace  applications,  v;ith  a  particular 
orientation  to  the  design-management,  decision-making  process.   Most  of 
the  papers  cited  treat  specialized  aerospace  structures  and  applica- 
tions; however,  the  author  does  try  to  generalize  by  introducing  a 
design  index  D,  a  material  efficiency  parameter  M,  and  a  structural 
efficiency  parameter  S.   After  defining  the  expressions  for  M  and  S 
corresponding  to  several  structural  elements,  design  charts  are  pre- 
sented wh.ich  show  regions  of  possible  application  for  various  materials. 
Unfortunately  the  design  charts  do  not  satisfy  expectations  aroused  by 
the  introduction  of  the  three  general  parameters. 

PvOzvany  (1966)  presents  a  similar  paper  pertaining  to  structures 
in  civil  engineering.   This  paper  is  less  comprehensive  and  more  ori- 
ented towards  specific  structural  applications  than  either  of  the  pre- 
ceding survey  papers.   The  author  postulates  several  "interrelated 
quantities  (or  parameters)"  which  could  perhaps  be  used  to  generalize 


structural  optimization  into  a  more  rational  methodology.   However, 
these  quantities — loading  (L) ,  material  (M) ,  geometry  (G) ,  initial 
behavior  (IB),  and  design  behavior  (SB) — are  only  applied  to  an  abstract 
discussion  of  concepts. 

Barnett  (1966)  has  a  readable  short  survey  of  the  field  that 
dwells  more  upon  theoretical  aspects.   He  postulates  a  general  problem 
in  which  the  cost  is  minimized  subject  to  a  system  in  equilibrium  with 
its  loads,  while  "behavior  constraints  on  strength,  stiffness,  and 
stability"  are  satisfied.   Uniform  strength  design  introduces  the  dis- 
cussion of  optimal  trusses;  virtual  work  theorems  that  were  derived 
originally  for  trusses  are  then  applied  to  simultaneous  plastic  collapse 
problems.   Following  a  brief  discussion  of  the  plastic  collapse  load 
bounding  theorems,  there  is  a  short  treatment  of  elastic  stability 
problems  and  material  merit  indices.   Barnett' s  stiffness  design  example 
exhibits  several  important  features  worth  noting.   Specifically,  the 
example  is  to  minimize  the  weight  of  a  beam  subject  to  some  given  load 
where  the  deflection  at  a  certain  point  is  specified.   Virtual  work  is 
used  to  handle  the  subsidiary  deflection  condition.   Necessary  condi- 
tions are  obtained  from  the  calculus  of  variations,  but  m.ore  signif- 
icantly, the  Schwarz  inequality  provides  a  sufficient  condition  for 
global  optimality.   Barnett  concludes  his  survey  with  a  section  stating 
that  multiload  designs  satisfying  "all  three  behavior  criteria"  are  more 
easily  solved  in  "design  space"  by  mathematical  programming  techniques 
described  by  Schmidt  (1966) . 

As  a  sequel  to  the  comprehensive  survey  (1638  to  1962)  by 
Wasiutynski  and  Brandt,  Sheu  and  Prager  (1968b)  present  a  complete  review 


of  developments  from  1962  to  1968.   This  paper  contains  three  major 
sections:   general  background,  methodology,  and  specific  problems. 
In  the  first  section  they  state  that  the  "well-posed  problem  of  optimal 
structural  design"  requires  specification  of  the 

(1)  purpose  of  the  structure  (load  and  environment), 

(2)  geometric  design  constraints  (limits  to  design 
parameters) , 

(3)  behavioral  design  constraints  (limits  to  the  "state"  of 
the  structure) , 

(4)  design  objective  (cost  functional) . 

The  methodology  section  is  not  noteworthy,  but  the  third  section  lists 
what  is  in  their  opinion  the  specific  problem  types:   static  compliance, 
dynamic  compliance,  buckling  load,  plastic  collapse  load,  multipurpose/ 
multiconstraint  structures,  optimal  laj'out  (e.g.,  trusses),  reinforced/ 
prestrcssed  structures,  and  from  the  background  section,  probabilistic 
problems.   Their  concluding  remarks  succinctly  summarize  the  paramount 
difficulty  of  the  subject:   realistic  problems  are  too  complicated  for 
precise  analytical  treatment.   \vhile  progress  is  being  made  in  analyt- 
ical treatment  of  simple  structures,  the  authors  opine  that  realistic 
problems  require  mathematical  programming  techniques.   Hov/ever,  they  do 
feel  that  analytical  treatment  is  desirable  to  provide  "a  deeper  insight 
into  the  analytical  nature  of  optimality." 

A  related  survey  by  Wang  (1968)  on  distributed  parameter  systems, 
".  .  .  whose  dynamical  behaviors  are  describable  by  partial  differential 
equations,  integral  equations,  or  functional  differential  equations," 
consists  entirely  of  a  bibliography.   While  not  pertinent  to  the  disser- 
tation, it  is  mentioned  here  for  completeness. 


Prager  (1970)  provides  another  survey  which  is  not  comprehensive, 
nor  does  he  present  any  new  results  as  claimed.   For  example,  Prager 
and  Taylor  (1968)  used  the  principle  of  minimum  potential  energy 
and  the  assumption  that  stiffness  is  proportional  to  specific  mass  to 
prove  global  optimality.   It  could  hardly  be  called  a  new  development 
in  1970.   At  the  same  time,  the  author  does  present  an  excellent  example 
of  a  multipurpose  optimal  design  problem.   Prager  also  treats  "segment- 
wise  constant"  approximations  and  the  optimal  layout  of  trusses. 

The  final  survey  paper,  Troitskii  (1971),  is  an  unusual  review 
of  methods  in  the  calculus  of  variations.   VJhereas  some  of  the  earlier 
surveys  present  lengthy  lists  of  references  but  contain  little  method- 
ology or  theory,  this  survey  is  just  the  opposite.   Part  of  what  makes 
it  unique  is  that  the  author  believed  only  eight  articles  merited  cita- 
tion— all  of  them  by  Troitskii.   This  shortcoming  is  more  than  overcome 
by  a  thorough  classification  of  optimal  control  problems  in  the  calculus 
of  variations.   Troitskii  bases  the  classification  on  "certain  character- 
istics of  control  problems":   types  of  constraints,  properties  of  the 
governing  dynamical  equations  of  the  system,  type  of  cost  functional,  and 
possible  state  discontinuities.   From  these  four  criteria  he  postulates 
five  principal  classes  of  problems;   however,  it  is  the  criteria  that 
are  important  and  not  the  specific  problem  type. 

By  comparing  what  the  authors  of  the  aforementioned  surveys 
believe  to  be  the  important  types  of  problems,  it  is  readily  apparent 
that  there  is  little  agreement  on  which  characteristics  of  structural 
optimization  problems  are  significant. 


1. 1   Historical  Development: 
Optimal  Columns 

The  beginning  of  structural  optimization  is  generally  attributed 
to  Galileo's  studies  in  1638  of  the  bending  strength  of  beams.   Accord- 
ing to  Barnett  (1968),  Galileo  considered  a  constant-width  cantilever 
beam  under  a  tip  load  as  part  of  a  study  of  "solids  of  equal  resistance," 
In  requiring  the  maximum  stress  in  each  cross  section  to  be  constant 
throughout  the  beam,  the  height  must  be  a  parabolic  function  of  position 
along  the  beam.   While  this  appears  to  be  the  origin  of  the  field,  a  prob- 
lem that  received  more  attention  is  the  buckling  of  a  column. 

Using  the  newly  developed  calculus  of  variations,  in  1773  Lagrange 
attempted  to  apply  variational  techniques  to  the  problem  of  finding  that 
distribution  of  a  homogeneous  m.aterial  along  the  length  of  a  column  which 
maximizes  the  buckling  load.   Truesdell  (1968)  relates  that  through  an 
insufficient  mathematical  formulation  Lagrange  shotved  the  optimal  form 
to  be  a  circular  cylinder.   Clausen  (185])  provides  the  earliest  known 
solution  to  this  problem  for  the  simply  supported  case.   As  described  in 
Todhunter  and  Pearson  (1893,  pp.  325-329),  Clausen  minimized  the  volume 
of  the  column  with  the  differential  equation  for  buckled  deflection 
treated  as  a  subsidiary  condition.   Assuming  all  cross  sections  to  be 
similar,  after  several  variable  transformations  and  complicated  manip- 
ulation, he  obtained  an  implicit,  analytical  solution. 

The  next  development  was  Greenhill  (1881),  according  to  Keller 
and  Niordson  (1966).   Greenhill  determined  the  height  of  a  uniform 
prismatic  column,  beyond  v.'hich  the  column  buckled  due  solely  to  its  ovm 
weight.   Timoshenko  and  Gere  (1961)  reproduce  the  solution  in  v/hich  the 


deflection  is  expressed  as  the  integral  of  a  Bessel  function  of  the 
first  kind  (of  the  negative  one-third  order). 

Blasius  (1913)  introduces  his  paper  with  a  uniform  strength  and 
a  minimum  deflection  beam  problem.   For  a  given  load  and  amount  of  mate- 
rial, the  cross-sectional  area  distribution  is  determined  which  maxi- 
mizes the  buckling  load  of  a  circular  column.   The  solution  is  identical 
to  that  obtained  by  Clausen.   In  addition,  Blasius  also  obtained  solu- 
tions for  columns  having  rectangular  cross  sections  and  discussed  the 
effect  of  different  boundary  conditions  on  the  results. 

For  the  next  few  decades,  structural  optimization  appears  to 
have  been  directed  towards  applications  in  the  aircraft  industry,  where 
aircraft  structural  problems  and  results  are  presented  in  the  format  of 
a  design  handbook.   Feigen  (1952)  is  a  good  example  of  this,  consider- 
ing the  buckling  of  a  thin-wall  column.   Given  a  constant  load  and  v/all 
thickness,  he  required  the  variable  inside  diameter  to  be  chosen  such 
that  the  buckling  load  is  maximized.   Wall  thickness  is  selected  to  make 
local  buckling  and  Euler  buckling  occur  at  the  same  load.   Solid  tapered 
columns  having  blunt  ends  are  also  treated  for  assumed  stiffness  dis- 
tributions. 

Renewed  interest  was  aroused  by  Keller  (1960),  who  examined  the 
problem  from  the  point   of  view  of  the  theory  of  elasticity,  and  in 
choosing  the  cross-sectional  shape  to  give  the  maximum  stiffness. 
Neglecting  the  weight  of  the  column,  he  obtained  via  the  former  that 
twisting  the  column  does  not  affect  the  buckling.   Of  all  convex  cross 
sections,  the  equilateral  triangle  is  shown  to  have  the  largest  second 


moment  of  area  relative  to  a  centroidal  axis.   Hence,  from  the  defin- 
ition of  buckling  load,  the  "best"  cross-sectional  shape  is  the  equi- 
lateral triangle.   Keller  also  obtained  Clausen's  implicit,  analytical 
solution.   Subsequently,  Tadjbakhsh  and  Keller  (1962)  generalized  the 
problem  to  a  general  eigenvalue  problem  and  boundary  conditions  subject 
to  a  subsidiary  equality  constraint.   The  latter  corresponds  to  spec- 
ifying the  volume  (or  weight)  of  material  to  be  distributed  in  an  opti- 
mal manner.   Using  the  Holder  inequality  they  demonstrate  global  opti- 
mality  of  the  eigenvalue  for  the  hinged-hinged  column. 

Keller  and  Niordson  (1966)  examine  the  case  of  a  vertical  column 
fixed  at  the  base,  subject  to  a  vertical  load  at  the  tip  and  the  column's 
oiJTi  weight.   It  is  also  assumed  that  all  cross  sections  are  similar. 
Their  approach  is  to  state  the  problem  as  a  simultaneous,  dual  optimiza- 
tion of  the  Raylelgh  quotient.   The  eigenvalue  is  minimized  with  respect 
to  the  eigenfunction  and  maximized  \-jith   respect  to  cross-sectional  area 
distribution.   Specifying  the  volume  of  m.aterial  available  is  treated 
as  a  subsidiary  equality  constraint.   From  the  maximum  lov;est  eigenvalue 
the  maximum  height  at  buckling  is  determined.   Solutions  are  obtained  by 
an  iterative  technique  employed  with  integral  equations. 

In  a  brief  note,  Taylor  (1967),  suggests  that  energy  methods  may 
link  cptimum  column  problems  to  the  traditional  eigenvalue  problems  of 
mechanics.   Prager  and  Taylor  (1968)  classify  problems  in  optimal  struc- 
tural design  and  demonstrate  global  optimality  using  energy  principles. 
Unfortunately  their  assumption  of  thin-v/all  construction  limits  the 
results  to  structures  vzhere  the  stiffness  is  proportional  to  the  specific 


mass  density.   The  consequence  of  this  assumption  is  that  in  the  energy 
formulation  the  resulting  control  law  and  governing  equations  are 
decoupled,  and  hence  easily  solved.   Huang  and  Sheu  (1968)  apply  this 
same  thin-wall  assumption  to  the  problem  treated  by  Keller  and  Niordson. 
However,  the  former  seek  the  maximum  end  load  instead  of  the  maximum 
height.   The  authors  also  attempt  to  limit  the  maximum  allowable  stress 
and  obtain  solutions  by  a  finite-difference  method.   Further  discussion 
of  sandwich  (thin-wall  construction)  columns  is  given  by  Taylor  and  Liu 
(1968) .   Basically,  this  paper  is  an  elaboration  of  the  techniques 
described  by  Prager  and  Taylor  when  applied  to  columns.   Extensive 
results  are  provided  for  various  cases. 

Post-buckling  behavior  for  columns  subject  to  conservative 
loads  is  considered  by  Gajewski  and  Zyczkowski  (1970).   A  nonconcerva- 
tive  problem  is  treated  by  Plaut  (1971b).  The  first  of  these  two  papers 
is  lengthy  but  is  much  too  narrow  in  scope  to  be  particularly  useful. 
In  the  second  paper,  the  Ritz  method  is  applied  to  an  energy  functional, 
obtaining  the  "best"  form  of  the  assumed  approximation  to  the  optimal 
solution. 

1 . 2  Historical  Development: 
Optimal  Static  Beams 

That  beam  problems  played  a  role  in  the  early  developments  of 
structural  optimization  has  already  been  indicated  in  the  preceding  sec- 
tion.  No  attempt  is  made  in  what  follows  to  present  a  complete  history, 
but  merely  to  outline  the  type  of  problems  that  have  been  considered. 
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Opatowski  (1944)  has  an  outstanding  paper  that  deals  with 
cantilever  beams  of  uniform  strength.   Besides  providing  numerous  refer- 
ences to  earlier  studies,  the  author  treats  the  problem  with  impres- 
sive mathematical  rigor.   The  beam  is  considered  to  deflect  under  its 
own  weight  and  a  transverse  tip  load;  bending  deflection  is  described 
by  a  Volterra  integral  equation  which  is  solved  exactly  for  various 
assumed  classes  of  variable  cross-sectional  geometry.   This  paper  is 
representative  of  earlier  papers  in  that  it  contains  extensive  analysis 
and  analytical  results,  but  little  numerical  data.   Barnett's  work  (1961) 
and  its  sequal  (1963a)  apply  the  calculus  of  variations  to  more  real- 
istic I-beams.   One  problem  considered  is  maximizing  the  weight  sub- 
ject to  general,  unspecified  loads  for  a  specified  deflection  at  a  given 
point.   The  Schwarz  inequality  is  used  to  derive  a  sufficient  condition 
for  global  optimality.   Also  included  is  a  comparison  of  uniform  strength 
beams  to  the  minimum  deflection  beam  for  several  different  cases  of 
applied  load  and  geometry.   The  paper  is  concluded  with  various  minimum 
weight  design  examples  in  which  both  bending  and  shear  stiffness  are  con- 
sidered. 

Haug  and  Kirmser's  (1967)  work  is  one  of  the  most  comprehensive 
studies  of  minimum  weight  beam  problem  published,   Wiile  it  may  succeed  in 
handling  any  conceivable  problem  and  in  employing  the  most  realistic 
stress  constraints,  this  very  generality  requires  so  many  variables  and 
conditions  that  the  mathematics  is  complicated  almost  beyond  reason. 

Another  study  of  minimum  weight  beams  (Huang  and  Tang,  1969)   is 
important  for  several  reasons.   By  dividing  the  beam  into  many  segments 
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having  constant  properties,  and  determining  the  necessary  conditions 
that  must  be  satisfied  by  every  segment,  it  appears  that  the  authors 
are  using  the  same  methods  that  were  used  to  solve  the  Brachistochrone 
problem  in  the  seventeenth  century.    In  their  treatment,  multiple 
loads  must  produce  a  specified  deflection  at  a  given  point;  these 
"segmentwise  constant"  approximations  and  multiple  load  problems  have 
recently  received  more  attention.   Of  further  interest  in  this  paper  is 
the  derivation  of  the  multiple  load  optimality  condition  using 
Pontryagin's  Maximum  Principle. 

V/hile  no  new  results  or  techniques  are  contained  in  Citron's 
(1969,  pp.  154-165),  the  author  gives  a  very  readable  minimum  weight 
beam  example.   The  problem  is  simple,  described  in  detail  completely, 
and  provides  an  excellent  example  of  how  control  theory  is  applied  to 
an  optimal  problem  in  structural  design.   An  analysis  of  intermediate 
beams  which  may  form  plastic  hinges  is  provided  by  Gjelsvik  (1971). 
It  shows  that  if  hinges  are  placed  at  all  points  of  the  beam  where  the 
bending  moment  is  zero,  this  makes  the  plastic  or  elastic  minimum 
weight  beam  statically  determinate.   Both  the  elastic  and  plastic  beam 
designs  are  shown  to  be  fully  stressed,  i.e.,  are  uniform  strength  beams, 

Application  of  generalized  vector  space  techniques  is  character- 
istic of  more  recent  papers.   Bhargava  and  Duf f in's  (1973)  is  such  a  study, 
It  treats  the  maximum  strength  of  a  cantilever  beam  on  an  elastic 


Since  these  early  papers  are  not  readily  available,  this  statement 
is  based  upon  descriptions  of  them  given  by  historical  references 
cited  in  Appendix  A. 
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foundation  subject  to  an  upper  bound  on  weight.   Although  involving 
more  advanced  mathematics  than  normally  required  by  variational  tech- 
niques, vector  space  methods  may  also  provide  more  powerful  analyt- 
ical tools. 

1. 3  Historical  Development: 
Optimal  Dynamical  Beams 

The  title  of  this  section  is  a  misnomer.   On  the  basis  of  the 
related  literature,  a  more  realistic  terminology  would  be  "Optimal 
Quasistatic  Beams."   Papers  dealing  with  the  optimization  of  djmamical 
beams  ultimately  contain  some  assumption  or  given  condition  that  effec- 
tively transforms  the  problem  to  an  equivalent  static  case.   Simple 
harmonic  motion  is  frequently  assumed  to  remove  time  dependence  from  the 
governing  equation.   For  example,  Barnett  (1963b)  minimizes  the  tip  of 
deflection  of  a  cantilever  beam  accelerating  uniformly  upwards,  with  the 
total  weight  of  the  beam  specified.   The  optimal  solution  for  cross- 
sectional  area  distribution  is  specified  by  a  nonlinear  integral  equa- 
tion solved  by  successive  approximations.   Dynamics  enters  the  problem 
only  as  a  time  invariant  inertia  load  vzhich  converts  the  problem  to  a 
static  beam  subject  to  a  body  force. 

Niordson  (1964)  finds  the  tapering  of  a  simply  supported  beam 
of  given  volume  and  length  which  maximizes  the  fundamental  frequency  of 
free  vibration.   Assuming  that  all  cross  sections  are  similar,  Niordson 
expresses  the  desired  frequency  as  the  Rayleigh  quotient.   This  is 
obtained  by  the  equation  for  spatial  dependence  associated  with  the 
usual  separation  of  variables  technique.   In  this  case  it  is  assumed 
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that  deflection,  shear,  and  rotational  inertia  effects  are  small. 
Solution  of  the  conditions  for  optimality  were  obtained  by  successive 
approximations.   This  approach  results  in  a  problem  identical  in  form 
to  the  eigenvalue  problem  associated  with  maximizing  the  Euler  buckling 
load  of  a  column. 

A  very  specialized  problem  is  treated  in  Brach  (1968).   Cross- 
sectional  area  and  stiffness  are  proportional,  and  have  both  upper  and 
lower'  bounds;  material  properties  and  length  of  the  beam  are  specified 
constants.   The  object  is  to  make  the  fundamental  frequency  of  free 
vibration  stationary,  without  any  constraints  on  the  weight.   Instead 
of  using  the  Rayleigh  quotient,  Brach  uses  the  total  potential  energy. 
His  solution  method  is  ad  hoc  and  not  generally  useful.   A  more  useful 
approach  is  described  by  Icerman  (1969)  for  structures  subject  to  per- 
iodic loads.   Necessary  and  sufficient  conditions  for  minimum  weight 
are  obtained  from  the  principle  of  minimum  potential  energy.   Amplitude 
and  frequency  of  the  applied  load  is  specified  as  well  as  "dynamic 
response,"  defined  as  the  potential  energy  associated  with  the  load 
amplitude  displaced  a  distance  equal  to  the  displacement  amplitude  at 
the  point  of  application.   It  is  also  required  that  the  load's  frequency 
be  less  than  the  fundamental  natural  frequency.   Subject  to  these  con- 
straints, the  structure's  weight  is  to  be  minimized.   Trusses  and 
segmentwise  constant  approximations  are  also  treated.   Once  again  a 
dynamical  problem  is  effectively  transformed  into  a  quasistatlc  problem. 

Realistic  treatment  of  the  dynamic  response  of  beams  subject  to 
time  dependent  loads  is  given  in  Plant  (1970).   An  upper  bound  on  the 
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response  of  the  beam  is  obtained  from  the  "largest  possible  displace- 
ment of  the  beam  under  a  static  concentrated  unit  load."   This  inequal- 
ity is  determined  from  the  time  derivative  of  the  total  energy  of  the 
beam  and  the  Schwarz  inequality.   Plaut  minimizes  the  upper  bound  on 
the  response  for  specified  total  weight  and  relationship  between  specific 
mass  and  stiffness.   Despite  the  consideration  of  a  truly  dynamical 
problem,  this  approach  has  two  vreaknesses  indicated  by  the  author. 
First,  the  upper  bound  need  not  be  close  to  the  exact  answer;  secondly, 
there  is  no  demonstration  that  minimizing  the  upper  bound  also  minimizes 
the  actual  response  of  the  beam  subject  to  dynamic  loading. 

Another  paper  worthy  of  mention  is  Brach  and  Walters  (1970) . 
They  maximize  the  fundamental  natural  frequency  which  is  expressed  as 
a  Rayleigh  quotient  that  includes  the  effect  of  shear.   Standard  varia- 
tional methods  are  employed  to  derive  necessary  conditions,  but  no 
solutions  nor  examples  are  given.   The  authors  do,  however,  suggest 
using  the  method  of  quasilinearization.   This  paper  is  another  example 
of  a  quasistatic,  dynamical  beam  problem. 

1 . 4   Scope  of  the  Dissertation 

Tills  dissertation  is  primarily  concerned  V7ith  the  application 
of  Pontryagin's  Maximum  Principle  to  problems  in  structural  optimization. 
Only  elastic  materials  are  considered;  however,  various  types  of  con- 
straints are  treated. 

The  theoretical  development  of  the  dissertation  pertains  to 
problems  described  by  an  ordinary  differential  equation  but  is  based 
upon  a  numerical  technique  normally  used  with  systems  described  by 
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a  finite  number  of  discrete  quantities.   For  this  reason,  examples  are 
included  for  both  types  of  systems. 

No  new  solution  techniques  have  been  developed  for  the  nonlinear 
two  point  boundary  value  problems  which  characteristically  arise  in 
optimization  problems  for  continuous  systems.   Solutions  are  obtained 
by  a  standard  quasilinearization  method.   However,  a  modified  "feasible 
direction"  numerical  algorithm  for  use  with  discrete  systems  is 
described  and  an  example  included  to  demonstrate  its  operation.   This 
serves  to  illustrate  the  application  of  the  theory  on  which  the  algo- 
rithm is  based  to  the  theoretical  development  associated  with  contin- 
uous systems.   Furthermore,  it  provides  a  comparison  between  the 
solution  to  a  problem  described  as  a  continuous  system,  and  alternately 
by  a  discrete  approximation. 

Additionally,  it  is  the  intent  of  this  dissertation  to  replace 
some  of  the  confusion  in  classification  of  problem  types  contained  in 
the  various  survey  papers   with  an  organization  based  upon  mathematical 
attributes.   The  result  is  a  logical  approach  to  the  formulation  of 
optimization  problems  for  elastic  structures. 


CHAPTER  IT 

GENERAL  PROBLEMS  AND  METHODS 
IN  STRUCTURAL  OPTIMIZATION 


2.0  Introduction 

The  first  chapter  presents  a  broad  view  of  structural  optimiza- 
tion and  the  historical  development  of  two  general  types  of  problems 
that  are  used  as  examples  in  later  chapters.   Before  the  mathematical 
theory  is  developed  in  the  next  chapters,  general  problems  and  methods 
in  structural  optimization  itself  are  briefly  outlined  in  this  chapter. 
The  classification  of  problem  types  vis-^-vis  mathematical  attributes 
is  discussed  first.   This  is  follo'^ed  by  short  descriptions  of  the 
major  methods  of  structural  optimization  for  both  continuous  and  dis- 
crete systems.   Appropriate  references  are  cited  for  each  section. 

2. 1  Problem  Classification  Criteria 

Perhaps  the  major  source  of  difficulty  in  classifying  struc- 
tural optimization  problems  lies  in  the  translation  from  the  physics 
involved  to  a  mathematical  representation.   A  single  physical  concept 
when  transformed  to  mathematics  may  become  more  than  one  mathematical 
attribute.   For  example,  consider  the  class  of  conservative  problems. 
Since  energy  is  conserved  this  immediately  prohibits  dissipative  mate- 
rials, time  varient  constraints,  and  nonholonomic  constants.   More 
important,  consider  the  following  statements  from  Lanczos  (1962, 
p.  226): 
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.  .  .  all  the  equations  of  mathematical  physics  which 
do  not  involve  any  energy  losses  are  deducible  from 
a  "principle  of  least  action,"  that  is  the  principle 
of  making  a  certain  scalar  quantity  a  minimum  or  a 
maximum  .  .  .  all  the  differential  equations  which  are 
self-adjoint,  are  deducible  from  a  minimum-maximum 
principle  ,  .  .  and  vice  versa. 

However,  it  is  shown  in  Chapter  V  that  to  be  self-adjoint  systems  places 
requirements  on  both  the  differential  equation  and  the  boundary  condi- 
tion.  Thus,  the  single  physics  attribute  of  being  a  conservative  sys- 
tem is  described  mathematically  by  expressions  involving  the  cost  func- 
tional, the  governing  system  of  differential  equations,  boundary  condi- 
tions, and  constraints. 

As  a  result  of  this  lack  of  similarity  in  descriptions,  a  choice 
must  be  m.ade  as  to  which  realm  will  be  used  for  classification  of  prob- 
lems.  Since  all  problems  are  ultimately  transformed  to  mathematics, 
mathematical  characteristics  are  selected  as  the  criteria.   On  the  basis 
of  survey  paper  contents  and  the  many  related  papers,  it  is  felt  that 
the  proper  (not  necessarily  the  best,  nor  all  inclusive)  characteristics 
for  classification  of  problem  types  are: 
(i)   cost  functional 

(ii)   system  equation  and  boundary  conditions 
(iii)   control  constraints 

(iv)   behavioral  constraints  (state  and/or  control  constraints) 
Subsequent  discussion  is  in  terms  of  these  four  characteristics. 
Exceptions  to  these  descriptors  are  readily  acknowledged,  e.g.,  whether 
the  problem  is  deterministic  or  probabilistic.   An  example  of  the  latter 
may  be  seen  in  Moses  and  Kinser  (1967) .   These  exceptions  do  not  serve 
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as  negating  counterexamples  but  instead  indicate  the  requirement  for 
additional  descriptors  and  verify  the  difficulty  of  the  task,  suggest- 
ing the  need  for  further  comprehensive  study. 

2.1.0  Problem  Classification  Guidelines 

In  the  following  sections  each  criterion  is  briefly  discussed. 
Some  of  the  various  characteristics  of  each  are  mentioned,  and  where 
appropriate  references  exist  the  citation  is  given.   It  must  again  be 
emphasized  that  the  following  is  not  all  inclusive;  it  is  an  attempt  to 
categorize  the  types  of  problems  existing  in  the  literature  according 
to  the  four  mathematical  descriptors  postulated.   Moreover,  the  descrip- 
tors are  not  discussed  in  the  order  given  but  instead  are  treated  in 
the  order  normally  encountered  during  a  problem  solution. 

2.1.1  Governing  Equations  of  the  System 

The  imm.ediate  question  to  be  answered  is  whether  the  structural 
system  is  described  by  a  set  of  continuous  functions  or  a  set  of  dis- 
crete coiistants.   Bending  deflection  of  simple  structural  elements  is 
an  example  of  the  former;  design  of  trusses  is  a  good  example  of  the 
latter.   In  general,  variational  techniques  are  employed  with  contin- 
uous systems  while  mathematical  programming  techniques  are  most  fre- 
quently applied  to  the  discrete  systems.   However,  variational  methods 
can  be  applied  to  the  approximation  of  continuous  systems  by  discrete 
elements.   This  usually  takes  the  form  of  either  a  finite  element  or 
"segment^^?ise  constant"  approximation  of  the  continuous  structure. 
References  treating  the  different  systems  described  above  are  included 
in  the  section  on  methods. 
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2.1.2   Constraints 

Two  of  the  descriptors  are  postulated  to  be  control  constraints 
and  behavioral  constraints.   A  further  consideration  is  whether  the 
constraint  is  defined  by  an  equality  relationship  or  by  an  inequality 
expression.   Equality  constraints  are  handled  by  a  long-known  technique 
entitled  "isoperimetric  constraints."   Valentine's  (1937)  work  is  known 
to  contain  the  initial  development  of  a  technique  for  converting  inequal- 
ity constraints  to  equality  constraints.  The  introduction  of  slack 
variables  increases  the  number  of  variables  in  the  problem  to  be  solved 
but  at  the  same  time  permits  all  of  the  isoperim.etric  techniques  to  be 
used.   A  detailed  application  of  this  approach  is  presented  in  Appendix  B. 
It  should  also  be  noted  that  isoperimetric  constraints  are  sometimes 
referred  to  as  "accessory  or  subsidiary  conditions." 

Most  real  structural  optimization  problems  possess  an  isoperi- 
metric constraint  as  well  as  inequality  constraints  dependent  upon  the 
control  u.  and/or  the  state  x  of  the  structural  system.    Typicall)'  the 
control  constraints  are  the  result  of  geometrical  limitations  or 
restrictions  to  the  types  of  available  materials.   Behavioral  constraints 
are  related  to  the  state  of  the  system  and  may  depend  solely  upon  the 
state  (of  deform.ation) ,  or  in  the  case  of  most  stress  constraints, 
jointly  upon  the  state  and  control  variables.   With  this  distinction  the 
constraints  may  be  classified  vis-^-vis  the  two  criteria  and  optimal 
control  characteristics  as  follows: 


By  convention  all  vectors  are  column  vectors  unless  indicated  otherwise. 
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(i)   unconstrained 

(ii)      <i>(u)        ^   0       control  constraint 

(iii)      (t)(x,u)  <  0-\ 

>     behavioral  constraints 
(iv)      (})(x)    <  Oj 

All  of  this  discussion  pertains  to  both  continuous  and  discrete  systems. 

No  references  for  unconstrained  optimization  problems  are  given. 
It  may  be  that  in  some  problems  the  unconstrained  structural  optimiza- 
tion solutions  have  either  infinite  stiffness  and  finite  weight,  or 
finite  stiffness  and  zero  weight.   A  discussion  of  this  can  be  found  in 
Salinas  (1968,  pp.  23-26). 

Investigation  of  control  constraints  led  to  the  development  of 
the  maximum  principle.   Although  reserving  discussion  of  the  method  for 
a  later  section,  the  classical  reference  detailing  the  derivation  of 
the  principle  is  given  here  for  completeness.   Rozonoer  (1959)  treats  con- 
trol constraints  but  only  as  related  to  the  development  of  the  maximum 
principle. 

Most  of  the  literature  concerns  either  bounded  control  problems 
or  a  more  general  form  of  constraint  which  can  be  classified  as  a  behav- 
ioral constraint.   The  latter  is  a  mixed  constraint  which  depends  upon 
both  the  state  and  control  variables.   Breakwell  (1959)  is  a  lucid  paper 
dealing  with  this  type  of  constraint.   References  which  treat  state 
and/or  mixed  constraints  are  Bryson  et  al.  (1963),  and  Speyer  and 
Bryson  (1968).   Constraints  which  depend  upon  only  the  state  are  not 
treated  in  this  dissertation;  an  example  of  such  a  constraint  is  to 
determine  the  optimal  solution  for  some  problem  subject  to  an  upper 
bound  on  deflection  of  the  structure  at  any  point. 
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2.1.3   Cost  Functionals 

There  are  two  basic  types  of  cost  functionals  that  occur  in  the 
field  of  structural  optimization.  They  were  first  identified  by  Prager 
(1969)  but  not  for  the  proper  reasons.   Using  Prager's  notation,  they  are 


J  =  Min  /   F(i|;)  dV 
V 

/  GOl>)    dV 
V 


J^  =  Min 


^  j   H(^)  dV 

V 

where  F,  G,  and  H  are  sca].ar  functionals  of   ^j.   The  latter  functional 
represents  a  Rayleigh  quotient  associated  with  an  eigenvalue  problem. 
It  can  be  reduced  to  the  first  type  of  functional  shown  above  by  choos- 
ing a  normalization  of  the  eigenfunction  such  that  the  numerator 
equals  unity  for  all  admissible  variations  (see  section  5.3).   This 
normalization  is  thereafter  treated  as  a  subsidiarj'  constraint. 

VJliat  actually  distinguishes  the  second   functional  from  the 
first  is  not  that  the  functional  is  a  quotient,  but  that  the  extrem- 
Ization  of  an  eigenvalue  requires  a  dual  extremization  (see  section  5.3) 
In  terms  of  a  state  x  and  control  _u,  the  fundamental  eigenvalue  is 
given  by  minimization  of  the  Rayleigh  quotient  v/ith  respect  to  the 
eigeiif unction  X'    °^ 


J„  =  Min 


/  G(x,u)  dV 
V 


^    X   I  H(x,u)  dV 
V 

v/here  u  represents  some  specified  design  parameter.   If  the  desired 

result  is  to  maximize  the  cost  J  with  respect  to  all  admissible  _u,  it 

is  observed  that  a  second  extremization  is  required;  for  example,  see 
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Keller  and  Niordson  (1966).   Thus,  a  more  appropriate  manner  of  classi- 
fying cost  functionals  is  on  the  basis  of  whether  the  problem  statement 
implies  a  single  or  a  double  extremization.   Hence  the  two  basic  types 
of  cost  functionals  encountered  in  structural  optimization  are 


J  =  Min  /   F(x,u)  dV 
u   V 


J   =  Max  Min: 


/   G(x,u)  dV 

V 


u    X  /   H(x,u)  dV 
V 

where  x   must  satisfy  an  equilibrium  condition  of  the  state,  and  _u  is 

subject  to  some  admissibility  requirements. 

There  is  a  special  case  related  to  these  two  in  which  the  weight 

is  to  be  minimized  for  a  specified  eigenvalue.   This  problem  is 

treated  in  Icerman  (1969)  i%?ith  a  mathematical  discussion  of  such  a 

variational  problem  presented  in  Irving  and  Mullineux  (1959,  p.  394). 

In  terms  of  the  two  cost  functionals,  the  special  case  is 

J  =  Min  /   {G(x,u)  -  J^H(x,u)}  dV 
u   V  ^ 

where  J   is  a  specified  constant.   This  approach  is  frequently  employed 
in  eigenvalue  problems  to  avoid  the  inherent  difficulties  associated 
with  the  dual  extremization  problem. 

There  have  also  been  many  papers  published  that  consider 
"multi-purpose  structures,"  e.g.,  Prager  (1969),  Martin  (1970),  and 
Prager  and  Shield  (1968).   The  cost  function  for  such  problems  is 
defined  as 
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J  =   E   a .  J .  (2i>H) 
i=l  ^  ^ 

where  a.  are  positive  constants,  serving  as  weighting  parameters. 
l-Jhile  perhaps  demonstrating  much  potential,  no  significant  results 
obtained  with  this  approach  have  so  far  been  published.   I-Jliat  problems 
have  been  solved  are  too  simple;  indeed  the  authors  indicate  the  need 
for  using  a  discrete  approximation  and  mathematical  programming  tech- 
niques in  realistic  applications. 

A  subject  closely  related  to  "multi-purpose  structures"  is 
that  of  multiple  constraints.   It  is  mentioned  here  only  because  most 
papers  on  the  latter  also  include  the  former — see  Martin  (1971) .   The 
idea  of  multiple  constraints  is  not  new;  in  both  variational  and  math- 
ematical programming  fields  there  exist  standard  techniques  for  han- 
dling multiple  constraints. 

A  recent  Russian  paper  (Salukvadze,  1971)   suggests  an  alter- 
nate to  the  "multi-purpose"  cost  functional.   Instead  of  treating  a 
vector  functional  that  requires  the  choosing  of  weighting  coefficients, 
it  is  suggested  that  the  several  functionals  be  combined  into  one. 
Given  a  system  and  vector  cost  functionals 

J.[u.]  =  J.(_x,u,t)      i  =  l,...,k 

Let  u„„   denote  the  optimal  solution  for  which  J.  assumes  the  optimal 
—OPT  ^  1 

value  on  the  trajectory  of  the  system.   For  each  of  the  J.  there  is  a 

different  u, „  .   These  k  values  J.  can  be  thought  of  as  components  of 
-HJ?T  1 
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a  vector  r  where 


t 


r   =  (J^ 


iL 


(1) 


■  (k)f 


For  any  arbitrary  u   the  result 


r[u]  =  {J^  [u]  ...  J,  [u]} 
1  —       k  — 


.T 


is  just  some  vector  functional. 

Vector  r^  represents  a  constant  point  in  the  space  of 

(J  ,...,J,  ).   Since  no  choice  of  u  can  optimize  all  of  the  J.  simul- 
1      k  —  1 

taneously,  that  is,  to  attain  the  point  _r   in  J. -space,  the  best  alter- 
native  is  to  minimize  the  distance  between  ^[u]    and  r;  .   That  distance 
is  defined  by  the  Euclidean  norm.   To  avoid  the  question  of  incon- 
sistent dimensions  the  functionals  are  reduced  to  dimensionless  form. 
Thus, 

2 


k 

J[u]  =  T. 
i=l 


J.  [u]  -  J. 
1  —     1 


(i) 

L-^PTJ 


J. 
1 


-^PT 


and  u^„„  is  that  function  u  which  minimizes  the  functional  J[u 
—OPT  —  — ' 

Mathematically  speaking. 


u^p,^  =  ARGMIN  {J[u]} 


This  type  of  vector  cost  functional  is  much  more  appealing 
than  the  type  treated  in  the  papers  on  multi-purpose  structures. 
It  also  suggests  an  entirely  new  field  of  study:   the  more  realistic 


t  T 

Superscript  "T",  e.g.,  _u  ,  denotes  the  transpose  of  the  vector  u. 
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choice  of  cost  functionals.   The  mathematics  of  a  problem  seldom 
accommodates  financial  considerations.   For  example,  the  design  which 
requires  the  least  material  may  reduce  the  cost  of  materials  at  an 
overxizhelming  expense  in  manufacturing  or  fabrication.   \<7hen  aesthetic 
appeal  and  environmental  impact  are  included — as  must  be  done  in  any 
real,  commercial  application — the  selection  of  an  appropriate  cost 
functional  is  an  almost  insurmountable  task.   However,  a  simple  exten- 
sion of  Salukvadze's  composite  cost  functional  may  reduce  the  difficul- 
ties to  operations  research  considerations. 

In  problems  where  it  is  desired  to  optimize  simultaneously 
several  different  functionals,  not  all  having  the  same  dimensions,  the 
concept  of  a  generalized  inner  product  may  prove  useful.   It  is  defined 
in  terms  of  a  metric  operator  A;  for  some  general  vector  ^ 
II   112  A  /    N      T 

and  symbol  "="  means  "is  defined  by."  With  reference  to  the  vector 
cost  functional,  A  represents  a  set  of  scale  factors  which  converts 
all  of  the  separate  cost  functionals  to  a  common  dimension.   This  is 
where  the  operations  research  enters — relating  material  expense  to 
fabrication  to  sociological  considerations  and  so  forth — to  determine 
the  metric  A.   For  a  vector  e_  whose  elements  are  functionals, 

where  r[u]  and  r   are  defined  above,  the  composite  cost  functional  is 


J[u]  =  (r.,Ac) 
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The  weakness  in  this  method  is  that  an  optimal  solution  must  be 
obtained  for  each  individual  cost  functional  prior  to  attempting  a 
solution  to  the  composite  problem.   Additionally,  some  of  the  more 
abstract  cost  objectives  may  be  difficult  to  quantify  in  a  meaningful 
manner.   Despite  these  shortcomings  this  approach  does  suggest  inter- 
esting applications. 

2.2   Methods:   Continuous  Systems 

The  problems  characterized  by  mathematical  functions,  in  con- 
trast to  those  represented  by  a  set  of  discrete  constants,  are  normally 
treated  by  variational  techniques.   Many  books  on  this  subject  have 
been  published;  the  better  authors  include  Elsgolc  (1961),  Gelfand  and 
Fomin  (1963),  Dreyfus  (1965),  Hestenes  (1966),  Denn  (1969),  Luenberger 
(1969),  and  Bryson  and  Ho  (1969).   An  excellent  summary  paper  is  avail- 
able in  Berg  (1962) . 

To  see  how  these  techniques  are  applied,  three  papers  are 
recommended.   The  first  is  Blasius  (1913),  which  provides  sufficient 
detail  and  explanation  to  make  it  quite  worthwhile.   Although  it  does 
include  several  examples  that  involve  subsidiary  conditions,  no  inequal- 
ity constraints  are  treated.   An  example  that  includes  inequality  con- 
straints to  the  control  variable  is  coontained  in  Brach  (1968) .   A  much 
more  general  application  of  variational  principles  is  presented  in 
Oden  and  Reddy  (1974).   In  this  paper  a  dual-complementary  variational 
principle  is  developed  for  a  particular  class  of  problems.   It  is  shown 
that  the  canonical  equations  obtained  are  the  Euler-Lagrange  equations 
for  a  certain  functional. 
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2.2.0  Special  Variational  Methods 

Besides  the  ordinary  variational  method,  more  specialized 
techniques  have  been  developed  to  the  point  where  they  are  recognized 
as  independent  methods  in  their  own  right.   In  the  following  sections 
these  methods  are  identified  and  a  number  of  representative  references 
given. 

2.2.1  Energy  Methods 

The  oldest  of  these  methods  is  the  energy  method.   It  originated 
with  the  principle  of  minimum  potential  energy,  and  was  later  extended 
to  include  the  representation  of  eigenvalues  through  the  energy  func- 
tional.  A  good  discussion  of  the  former  is  available  in  Fung  (1965)  or 
Przemieniecki  (1968);  the  best  general  treatment  of  the  latter  is  avail- 
able in  either  Gould  (1957,  Chapter  4)  or  in  Mikhlin  and  Smolitskiy 
(1967,  Chapter  3). 

The  principle  of  minimum  potential  energy  is  frequently  used 
with  simple  problems  to  prove  that  a  necessary  condition  for  optimality 
is  also  sufficient.   Prager  and  Taylor  (1968)  contains  such  a  proof 
for  the  global,  maximum  stiffness  design  of  an  elastic  structure  of 
sandwich  construction;  two  papers  that  also  consider  this  problem  are 
Huang  (1968)  and  Taylor  (1969).   Specific  application  of  the  energy 
method  to  an  eigenvalue  problem  is  demonstrated  in  Taylor  and  Liu  (1968) . 
A  much  more  general  discussion  of  the  energy  method  is  provided  in 
Salinas  (1968) .   Further  extensions  of  the  method  are  presented  in 
Masur  (1970),  in  which  the  principle  of  minimum  complementary  energy 
is  applied  to  problems  of  the  optimium  stiffness  and  strength  of  elastic 
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structures.   In  these  problems  a   necessary  condition  for  optimality 
is  that  the  strain  energy  density  be  constant  throughout  the  structure. 
This  condition  is  also  sufficient  for  optimality  in  certain  classes  of 
structures  that  satisfy  a  specific  relationship  between  the  strain 
energy  density  and  design  variables. 

Many  of  the  energy  problems  belong  to  the  class  of  problems 
having  quadratic  cost  functionals.   The  significance  of  this  character- 
istic is  that  the  Euler-Lagrange  equations  derived  from  such  functionals 
are  linear. 

A  more  recent  energy  method  development  is  the  concept  of 
"mutual  potential  energy."   Mutual  potential  energy  techniques  resemble 
those  of  the  principle  of  minimum  potential  energy.   In  both  methods 
a  cost  functional  is  defined  over  the  entire  domain  occupied  by  the 
structure  and  optimized  with  respect  to  the  control  variable.   If  it 
is  desired  that  the  optimal  solution  be  required  to  have  a  specified 
deflection  at  a  certain  point,  this  condition  corresponds  to  a  sub- 
sidiary state  constraint  when  using  the  principle  of  minimum  potential 
energy.   The  mutual  potential  energy  method  incorporates  this  type  of 
localized  constraint  into  the  cost  functional  which  is  defined  over 
the  entire  domain  of  the  structure.   By  itself  this  alone  is  advanta- 
geous; however,  for  certain  types  of  problems  the  mutual  potential 
energy  method  also  provides  both  a  necessary  and  sufficient  condition 
for  global  optimality.   In  the  way  of  a  critical  comment,  either  the 
method  has  received  too  little  attention,  or  else  it  does  not  effi- 
ciently handle  problems  more  difficult  than  the  simple  examples  presented. 
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Four  papers  that  are  representative  of  the  literature  associated  with 
this  method  are  Shield  and  Prager  (1970),  Chern  (1971a,  1971b),  and 
Plaut  (1971c). 

Another  recent  development  is  that  the  class  of  problems  for 
which  energy  methods  is  applicable  has  been  expanded  to  include  cer- 
tain types  of  nonconservative  systems.   Together  with  the  mutual  poten- 
tial energy  concepts,  this  suggests  that  perhaps  the  classical  energy 
method  is  a  special  case  of  a  more  generalized  method.   If  a  technique 
can  be  developed  which  uses  the  adjoint  variables  to  transform  a  gen- 
eral nonconservative  system  into  an  equivalent  self-adjoint  form,  the 
method  might  be  deduced.   Some  papers  pertaining   to  the  subject  are 
Prasad  and  Herrmann  (1959),  Wu  (1973),  and  Barston  (1974). 

2.2.2   Pontryagin's  Maximum  Principle 

There  are  many  textbooks  ^^7hich  derive,  explain,  and  give  examples 
for  the  maximum  principle.   The  original  (Pontryagin  et  al.,  1962) 
requires  a  knowledge  of  functional  analysis.   A  condensed  form  of  this 
same  material  is  available  in  Rozonoer  (1959).   Denn  (1969)  provides 
another  point  of  view  in  which  the  principle  is  derived  from  Green's 
functions.   In  this  manner,  the  sensitivity  to  variations  is  readily 
observed.   To  understand  Denn's  treatment  requires   only  a  knowledge  of 
the  solution  of  differential  equations. 

Shortly  after  Pontryagin's  book  was  published,  many  papers 
devoted  to  the  theoretical  aspects  of  the  maximum  principle  were  pub- 
lished.  Some  of  the  more  readable  ones  are  Kopp  (1962,  1963),  Roxin 
(1963),  and  Halkin  (1963).   Another  early  paper  (Breakwell,  1959) 
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appears  to  be  a  completely  independent  derivation  of  the  maximum  prin- 
ciple.  Although  quite  general  in  the  mathematical  sense,  the  examples 
presented  are  trajectory  optimization  problems  and  not  a  general  type 
of  mathematical  problem.   This  may  be  an  explanation  for  what  seems  to 
be  a  lack  of  recognition  for  a  significant  achievement. 

The  application  of  PMP  to  problems  in  structural  optimization 
is  relatively  recent.   When  the  method  is  used,  one  of  two  difficulties 
is  often  encountered.   The  first  is  errors  in  the  formulation  of  the 
optimal  control  problem;  the  second  is  that  once  a  well-posed,  non- 
linear two-point  boundary  value  problem  (TPBVP)  is  obtained,  it  is  dif- 
ficult to  solve.   An  example  of  the  first  is  provided  by  Dixon  (1968) — 
the  correction  was  given  in  Boykin  and  Sierakowski  (1972).   De  Silva 
(1972)  provides  a  clear  presentation  on  the  application  of  PMP  to  a 
specific  problem,  but  includes  no  data  because  a  solution  could  not  be 
obtained.   Despite  the  failure  to  determine  the  solution,  this  paper 
is  worthwhile  for  its  lucid  discussion  of  the  Pl-fP  application.   Another 
paper  that  gives  a  good  specific  application  of  PMP  is  Maday  (1973). 
Although  much  analysis  is  presented  very  little  is  said  regarding  the 
solution  techniques. 

All  of  the  above  references  are  applicable  only  to  systems 
that  are  described  by  ordinary  differential  equations,  in  contrast  with 
the  calculus  of  variations  which  also  handles  problems  described  by 
partial  differential  equations.   Since  many  of  the  problems  of  mathe- 
matical physics  involve  partial  differential  equations,  an  extension  of 
PMP  to  include  this  class  of  problems  is  the  next  logical  development. 
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Some  work  has  already  been  done,  for  example,  Barnes  (1971)  and 
Komkov  (1972).   A  survey  of  these  "distributed  parameter  systems" — 
see  Section  1.0 — is  presented  in  Wang  (1968). 

2.2.3  Method  of  Steepest  Ascent/Descent 

This  method  is  frequently  cited  in  the  literature  for  trajec- 
tory optimization,  and  occasionally  in  references  related  to  optimal 
structures.   l-Then  the  method  of  quasilinearization  converged  for  the 
dissertation  example  problems,  there  was  no  need  to  investigate  other 
methods  such  as  the  method  of  steepest  ascent.   Consequently,  little  is 
said  about  it.   According  to  the  references,  it  is  applied  in  a  straight- 
forx^.'ard  manner.   Furthermore,  the  example  problem  solutions  presented 
seem  to  be  real  problems  and  not  academically  simple.   The  following 
four  papers  treat  the  method  in  general  with  trajectory  optimization 
applications:   Bryson  and  Denham  (1962,  1964),  Bryson  et  al.  (1963),  and 
Hanson  (1968).   In  Haug  et  al.  (1969),  the  method  of  steepest  ascent  is 
derived  in  detail,  completely  discussed,  and  compared  to  the  maximi.mi 
principle.   Several  structural  optimization  problems  are  then  solved 
by  the  method  of  steepest  ascent.   Although  no  exciting  results  are 
obtained,  use  of  the  method  is  clearly  illustrated  by  the  applications 
to  realistic  structural  problems. 

2.2.4  Transition  Matrix: 
Aeroelasticity  Problems 

For  the  past  few  years  a  group  at  Stanford  University  has  studied 
the  optimization  of  structures  subject  to  dynamic  or  aerodynamic  con- 
straints.  The  general  problem  of  their  interest  is  that  of  minimizing 
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the  weight  of  a  given  structure  for  specified  eigenvalue,  subject  to 
inequality  constraints  on  control. 

Three  types  of  solution  techniques  are  used  once  the  necessary 
conditions  for  minimum  weight  are  determined.   Exact  solutions  are 
obtained  for  most  of  the  problems  because  they  are  so  simple  that 
analytical  methods  are  applicable.   More  complicated  problems  are 
solved  by  a  "transition  matrix"  method  described  in  Bryson  and  Ho  (1969)  . 
On  the  basis  of  convergence  difficulties  reported  in  the  references, 
this  method  should  be  used  with  caution.   Results  have  been  obtained 
only  for  very  simple  problems.   However,  these  results  are  corroborated 
by  data  obtained  from  a  discrete  approximation  method.   Five  papers 
that  are  representative  of  this  work  are  Mcintosh  and  Eastep  (1968), 
Ashley  and  Mcintosh  (1968),  Mcintosh  et  al.  (1969),  Ashley  et  al .  (1970), 
and  Weisshaar  (1970) . 

2.2.5   Miscellaneous  Methods 

The  preceding  sections  have  briefly  outlined  the  methods  of 
structural  optimization  most  frequently  encountered  in  the  literature. 
Appropriate,  representative  references  have  also  been  given.   Not  all 
methods  are  listed;  while  some  are  omitted  for  not  being  generally  use- 
ful, othBrs  are  omitted  for  not  being  generally  used.   T\v70  examples  of 
the  latter  are  the  "modified  quasilinearization"  and  "sequential  gradient- 
restoration"  algorithms  described  in  Miele  et  al.  (1972)  and  in  Hennig 
and  Miele  (1972) .   At  some  later  time  these  methods  may  be  acknowledged 
as  major  methods  that  are  applicable  to  many  different  or  important  prob- 
lems, but  for  now  chey  are  mentioned  only  in  passing. 
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2 . 3  Methods:   Discrete  Systems 

Discrete  systems  are  described  by  a  set  of  discrete  constants 
instead  of  the  set  of  functions  associated  with  continuous  systems. 
The  classic  example  of  a  discrete  system  is  a  pin-connected  truss,  where 
the  state  x.  and  control  u.  are  the  stress  and  cross-sectional  area, 

X  1 

respectively,  for  each  member  i.   A  discrete  system  also  arises  in  the 
approximation  of  a  continuous  system. 

Several  references  that  present  a  good  discussion  of  general 
methods  applied  to  discrete  systems  are  available.   Most  of  these  exist 
in  the  form  of  an  edited  collection  of  papers  by  various  authors  on  the 
topics  of  their  acknowledged  expertise.   Four  such  publications  are 
Gellatly  (1970),  Gellatly  and  Berke  (1971),  Pope  and  Schmidt  (1971), 
and  Gallagher  and  Zienkiewicz  (1973) .   Another  report,  Melosh  and  Luik 
(1967),  provides  a  good  e:cposition  of  the  difficulties  associated  with 
the  analysis  portion  of  least  weight  structural  design.   It  also  con- 
tains a  brief  comparison  of  various  mathematical  programming  methods. 

McNeill  (1971)  is  the  last  reference  to  be  cited  in  the  section 
on  general  methods  for  the  optimization  of  discrete  systems.   Minimum 
weight  design  of  general  structures  is  treated  in  a  mathematically 
precise  formulation.   Legendre's  necessary  condition  is  combined  with 
the  concepts  of  convex  functions  and  sets  to  derive  the  necessary  and 
sufficient  conditions  for  global  optimality.   Fully  stressed  designs 
and  constraints  to  eigenvalues  are  also  discussed.   In  summary,  this 
paper  provides  a  good  exam.ple  of  the  general  mathematical  problem  that 
must  be  solved  in  the  optimization  of  discrete  systems. 
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l^ile  certain  variational  methods  may  be  applied  to  discrete 
systems,  the  most  frequently  used  technique  is  mathematical  programming. 
In  the  following  sections,  this  method  and  other  major  methods  are  dis- 
cussed and  representative  references  cited. 

2.3.0  Mathematical  Programming 

The  general  method  of  mathematical  programming  is  discussed 
in  section  3.2  of  the  dissertation,  and  the  solution  of  an  example 
problem  using  this  method  is  detailed  in  Chapter  VL.  In  the  literature 
related  to  this  subject,  a  very  readable  textbook  is  available — Fox 
(1971) .   This  book  complements  the  theory  with  numerous  discussions 
pertaining  to  numerical  techniques  and  methods  that  can  be  employed  to 
overcome  certain  difficulties  that  may  arise.   Although  it  does  contain 
flowcharts  of  several  algorithms,  there  are  few  specific  examples  given. 
For  a  discussion  of  the  general  theory,  two  alternatives  to  this  book 
exist  in  the  form  of  papers:   Schmidt  (1966,  J968).   The  first  is 
written  in  a  conversational  style,  contains  no  mathem.atics,  and  is 
intended  to  provide  only  a  general  description  of  the  subject.   The 
latter  paper  is  theoretical  in  content. 

An  excellent  application  to  a  realistic  problem  is  to  be  found 
in  Stroud  et  al.  (1971).   Ths  paper  contains  little  discussion  of  the 
method  itself,  but  does  demonstrate  an  application  that  allows  a  concep- 
tual visualization  of  the  solution.   The  approach  is  to  assume  the 
solution  to  be  a  linear  combination  of  specified  functions,  and  to 
choose  the  weighting  coefficients  to  minimize  the  cost.   Mathematical 
programming  is  employed  to  determine  the  optimal  set  of  coefficients. 
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This  approach  resembles  Galerkin's  method,  and  though  not  mathemat- 
ically rigorous,  it  may  provide  a  useful  approximation  to  large, 
unwieldy  problems. 

2.3.1   Discrete  Solution  Approximations 

In  the  previous  section  a  paper  is  cited  that  contains  an 
approximate  solution  obtained  by  Galerkin's  method.   The  use  of  the 
Galerkin  or  Rayleigh-Ritz  approximate  solution  techniques  is  suffi- 
ciently widespread  to  be  considered  a  general  method.   For  both  methods, 
the  solution  is  assumed  to  be  a  linear  combination  of  the  solution  to 
the  linear  part  of  the  governing  equation  and  a  set  of  prescribed  func- 
tions.  This  approximate  solution  does  not  satisfy  the  given  equation 
exactly  but  produces  some  residual  function.   A  cost  function  that 
depends  upon  the  residual  is  then  minimized  \jith  respect  to  the  unknox^n 
coefficients.   The  two  weighted  residual  methods  mentioned  above  have 
different  cost  functions,  but  the  methods  are  identical  for  linear 
equations — see  Cunningham  (1958,  p.  158). 

The  advantage  to  using  these  methods  is  that  after  assuming 
the  particular  form  of  solution,  the  problem  of  solving  for  the  weight- 
ing coefficients  may  be  much  simpler  than  the  original  problem.   In  the 
case  of  Stroud  et  al.  (1971),  the  coefficients  were  obtained  by  mathe- 
matical programming  techniques.   However,  the  weakness  of  the  method 
is  the  restricted  function  space  of  possible  solutions.   With  the  coef- 
ficents  obtained  by  these  methods  the  resulting  solution  is  the  best 
approximation  that  is  possible  from  the  set  of  solution  functions 
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prescribed.   There  is  no  guarantee  that  the  approximation  even  resembles 
the  true  solution. 

The  flutter  of  a  panel  is  solved  using  Galerkin's  method  in 
Plaut  (1971a) .   No  general  developments  are  presented  and  the  assumed 
solution  functions  are  trivially  simple.   However,  this  paper  does 
provide  an  application  of  the  method  to  obtain  an  approximate  solution 
to  a  very  difficult  optimization  problem  involving  the  stability  of 
a  nonconservative  system.   A  similar  problem  is  treated  in  a  more  theo- 
retical manner  in  Plaut  (1971b)  using  a  modified  Rayleigh-Ritz  method. 
"Segmentwise-constant"  control  functions  are  assumed  also;  this  partic- 
ular approximation  is  discussed  with  more  detail  in  the  following  sec- 
tion.  Additional  nonconservative  problems  are  treated  in  Leipholz 
(1972),  applying  Galerkin's  approximate  solution  to  the  energy  method. 
Several  simple  examples  are  included. 

Nonconservative  elastic  stability  problems  of  elastic  continua 
are  treated  in  Prasad  and  Herrmann  (1969)  using  adjoint  systems.   This 
approach  is  more  realistic  than  the  segmentwise-constant  control  assump- 
tion described  in  the  following  section.   Solutions  for  the  state  and 
adjoint  system  are  assumed,  such  that  approximation  process  resembles 
the  Rayleigh-Ritz  method.   However,  only  a  single  type  of  nonconserva- 
tive system  is  considered.   Extension  to  several  other  types  of  noncon- 
servative elastic  continua  problems  is  given  in  Dubey  (1970).   Varia- 
tional equations  corresponding  to  both  the  Galerkin  and  Rayleigh-Ritz 
methods  are  derived.   Furthermore,  the  condition  for  equivalence  of  the 
two  methods  is  shown  to  be  that  the  admissible  velocity  field  must 
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satisfy  a  natural  boundary  condition  over  that  portion  of  the  body's 
surface  where  tractions  are  prescribed. 

2.3.2   Segmentwise- Constant  Approximations 

The  definitive  characteristic  of  this  method  is  approximating 
the  structural  system  by  a  number  of  discrete  segments,  where  within 
each  segment  the  control  function  has  a  constant  value.   In  general, 
the  constant  value  of  the  control  differs  from  segment  to  segment.   For 
the  many  papers  on  this  method  that  have  been  published,  the  procedure 
is  the  same.   An  optimality  condition  (necessary  in  all  cases  but  also 
sufficient  in  some)  or  cost  functional  is  derived  for  the  continuous 
system.   After  defining  the  segmentwise-constant  approximation,  the 
condition  or  functional  is  reformulated  in  terms  of  the  discrete  param- 
eters.  Most  of  the  papers  use  so  few  elements  that  solving  for  the 
discrete  values  of  the  control  parameter  poses  no  difficulties.   Although 
this  method  does  simplify  the  mathem.atical  problem  to  be  solved,  the 
crudeness  of  the  approximation  is  not  appealing.   Five  papers  which 
treat  a  variety  of  problems  using  this  approximation  are  cited  below. 

Miniraiim  weight  of  sandwich  structures  subject  to  static  loads 
is  discussed  in  Sheu  and  Prager  (1968a) .   In  Sheu  (1968)  the  same  type 
of  structure  is  considered.   It  differs  from  the  first  problem  by 
requiring  point  masses  to  be  supported  such  that  the  total  structure  has 
a  prescribed  fundamental  frequency  of  free  vibration.   Icerman  (1969) 
treats  the  problem  of  elastic  structures  subject  to  a  concentrated  load 
of  harmonically  varying  amplitude.   The  minimum  weight  design  is  obtained 
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subject  to  a  compliance  constraint  related  to  the  applied  load,  and 
which  is  effectively  a  boundary  condition  on  displacement  at  the  point 
of  application.   A  truss  problem  is  also  included. 

The  concept  of  a  compliance  constraint  is  pursued  further  in 
Chem  and  Prager  (1970).   The  minimum  weight  design  for  sandwich  con- 
struction beams  under  alternative  loads  is  found,  subject  to  this  type 
of  constraint.   The  paper  uses  up  to  eight  segments,  thereby  obtaining 
a  more  realistic  approximation  to  the  continuous  problem.   Minimum  weight 
design  of  elastic  structures  subject  to  body  forces  and  a  prescribed 
deflection  is  discussed  in  Chern  (1971a).   Tnis  investigation  is  no- 
table in  that  it  considers  applied  loads  that  are  functions  of  the  design 
functions. 

2.3.3   Complex  Structures  with 
Frequency  Constraints 

On  the  basis  of  useful  application,  perhaps  the  most  important 
class  of  discrete  structural  optimization  problems  is  the  minimum  weight 
design  of  complex  structures  subject  to  natural  frequency  constraints. 
Since  most  real  structures  are  built  with  many  structural  elements  of 
various  types,  and  are  not  realistically  described  by  any  single  type, 
this  approach  is  more  appropriate  from  the  aspect  of  modeling  the  struc- 
ture.  Furthermore,  many  structures  must  be  designed  to  avoid  certain 
natural  frequencies  because  of  resonance  or  self-induced  oscillations; 
this  situation  indicates  that  the  natural  frequency  constraint  is  also 
appropriate. 

Many  different  solution  schemes  have  been  developed  which  are 
usually  based  upon  general  mathematical  programming  techniques. 
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Typically,  a  design  is  iteratively  altered  to  minimize  the  weight  with 
a  subsequent  increase  in  frequency  until  a  constraint  is  violated. 
At  that  point  the  design  process  uses  an  iteration  which  simultaneously 
reduces  both  weight  and  frequency.   These  two  processes  are  repeated 
sequentially  until  no  further  weight  reduction  is  possible. 

Although  circumstances  may  require  the  use  of  many  elements, 
the  number  of  them  may  itself  be  a  critical  factor.   Some  of  the  schemes 
require  a  matrix  inversion  as  part  of  the  eigenvalue  problem  solution 
associated  with  the  frequency  constraint.   If  the  number  of  elements 
becomes  too  large,  the  size  of  the  matrix  to  be  inverted  likewise  be- 
comes excessively  large.   When  that  occurs  the  matrix  inversion  can 
require  excessive  amounts  of  computer  time.   Another  possible  difficulty 
is  that  the  inverse  matrix  itself  is  not  sufficiently  accurate,  such 
that  the  subsequent  calculations  are  not  acceptable.   However,  for 
structures  such  as  reinforced  shells  composed  of  different  types  of 
structural  elements,  this  method  may  be  the  most  applicable. 

Many  papers  have  been  published  pertaining  to  this  class  of 
structural  optimization  problem.   Because  the  method  is  inherently 
oriented  towards  applications,  the  references  are  cited  in  chronolog- 
ical order  without  additional  comments.   Interested  readers  are  referred 
to:   Turner  (1967),  Zarghamee  (1968),  Turner  (1969),  De  Silva  (1969), 
Rubin  (1970),  Fox  and  Kapoor  (1970),  McCart  et  al.  (1970),  and  Rudisill 
and  Bhatia  (1971). 
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2.3.4   Finite  Element  Approximations 

There  is  an  unfortunate  ambiguity  to  the  label  "finite  elements" 
that  occurs  because  these  words  are  used  to  describe  two  completely 
different  entities.   In  papers  cited  in  the  preceding  section  they  are 
used  to  indicate  the  discrete  structural  elements  of  finite  dimensions 
which  comprise  the  complex  structure.   The  analysis  of  such  systems  of 
structural  elements  has  been  accomplished  by  ordinary  matrix  methods 
during  the  last  three  decades.   However,  during  the  past  decade  another 
method  has  been  developed  and  named  "the  finite  element  method." 

In  this  method  a  continuum  is  divided  into  small,  finite  ele- 
ments over  which  a  particular  form  of  approximation  of  either  the  dis- 
placement and/or  force  is  assumed.   A  number  of  nodes  common  to  one  or 
more  element  is  prescribed;  continuity  is  required  to  exist  at  these 
nodes  but  not  necessarily  elsewhere.   An  equilibrium  equation  is  derived 
for  each  element,  and  then  all  of  the  individual  equations  are  combined 
into  a  single  equilibrium  equation  for  the  entire  system.   The  result- 
ing equation  is  a  linear  algebraic  equation  whose  unknowns  are  displace- 
ments and/or  forces  at  the  nodes.   Once  the  matrix  equation  is  inverted, 
the  nodal  displacements  and/or  forces  are  used  with  the  assumed  approx- 
imation form  to  describe  the  stats  of  the  structure  throughout  each  and 
every  elem?ni:,  and  hence  the  system.   Hereafter  this  method  is  referred 
to  as  tlie  "finite  element  method." 

The  most  frequent  application  of  the  finite  element  method  is 
to  problems  having  complicated  loads,  geometry,  and  response.   Generally 
speaking,  the  method  is  employed  wherever  the  physical  system  is  too 
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complex  to  be  described  adequately  by  a  single  differential  equation 
and  boundary  conditions.   For  a  complete  theoretical  development  of  the 
finite  element  method  and  numerous  examples,  see  Zienkiewicz  (1971). 

With  respect  to  structural  optimization  the  method  is  employed 
to  simplify  the  problem  to  be  solved.   Very  little  has  been  published 
on  this  subject,  but  the  papers  available  cover  a  wide  spectrum  of  tech- 
niques.  For  example,  Dupuis  (1971)  combines  the  finite  element  and 
segmentwise-constant  methods  as  applied  to  minimum  weight  beam  design. 
A  similar  application  to  column  buckling  is  contained  in  Simitses  et  al, 
(1973).   Another  paper,  Wu  (1973),  is  a  study  of  two  classical  noncon- 
servative  stability  problems.   Although  adapted  to  stability  consider- 
ations, this  presentation  is  the  best  exposition  available  in  the  open 
literature. 

In  Chapter  VI  a  minimum  deflection  beam  problem  is  solved  v/ith 
the  combined  methods  of  finite  elements  and  mathematical  programming. 

2 . 4   Closure 

In  the  preceding  sections  of  this  chapter,  general  problem 
types  and  methods  are  discussed.   Only  those  m.ethods  that  appear  to 
have  attained  some  standard  of  acceptance  are  presented.   It  must  be 
acknowledged  that  other  areas  of  important  study  exist  but  are  perhaps 
overlooked  as  not  being  pertinent  to  the  general  subject  area  of  the 
dissertation.   As  an  example,  Dorn  et  al.  (1964)  treats  the  optimal 
layout  of  trusses — an  important  subject  but  not  related  to  the  general 
problem  to  be  considered  in  this  dissertation.   In  addition  only 
elastic  structures  have  been  considered  although  there  are  numerous 
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publications  on  optimal  design  of  inelastic  structures.   References 
that  are  representative  of  this  subject  are:   Drucker  and  Shield 
(1957a,  1957b),  Hu  and  Shield  (1961),  Shield  (1963),  Prager  and  Shield 
(1967),  and  Mayeda  and  Prager  (1967). 

On  considering  the  various  references  mentioned  above  it  would 
appear  that  there  are  two  possible  pitfalls  in  structural  optimization 
that  should  be  avoided.   The  first  is  the  confusing  of  method  of  opti- 
mization with  the  solution  techniques  employed  to  obtain  a  solution  to 
the  resulting  TPBVP.   In  order  to  avoid  possible  errors  the  two  should 
be  dealt  with  Independently,  unless  it  is  clearly  advantageous  to 
relate  one  to  the  other.   Besides  this  it  must  be  recognized  that  any 
solution  obtained  is  "optimal"  only  with  respect  to  the  given  condi- 
tions of  the  particular  problem.   Any  change  in  the  problem  statement 
invalidates  the  applicability  of  that  solution.   The  change  may  lead  to 
a  more  desirable  solution,  but  the  original  solution  is  no  less  valid. 
Simitses  (1973)  is  an  example  where  this  situation  is  not  acknowledged. 
In  this  paper  the  thickness  of  a  thin  reinforced  circular  plate  of  spec- 
ified weight  and  diameter  is  determined  such  that  the  average  deflec- 
tion due  to  a  uniform  load  is  minimized.   An  earlier  paper  which  did 
not  include  stiffening  is  cited  with  the  implication  that  the  optimum 
solution  for  the  unstiffened  plate  is  not  correct.   The  point  made  above 
is  that  both  of  these  solutions  are  optimum  under  the  respective  condi- 
tions of  the  two  problems.   Neither  solution  is  more,  or  less,  valid 
than  the  other. 


CHAPTER  III 
THEORETICAL  DEVELOPMENT 

3.0   Introduction 

This  chapter  contains  the  development  of  two  distinct  methods 
used  in  the  theory  of  optimal  processes,  into  a  more  general  method. 
The  first  section  defines  precisely  the  problem  to  be  considered. 
This  includes  the  necessary  conditions  for  an  optimal  solution  given  by 
the  calculus  of  variations.   Several  mathematical  programming  techniques 
are  described  in  the  second  section  along  vzith  a  numerical  algorithm 
called  the  gradient  projection  method.   The  application  of  this  numer- 
ical method  to  the  solution  of  the  necessary  conditions  from  Pontryagin's 
Maximum  Principle  (Pr-lP)  is  detailed  in  Section  3.3.   Results  of  this 
approach  are  shown  to  be  consistent  with  the  necessary  conditions,  given 
in  Section  3.1;  these  results  provide  a  clarifying  insight  to  the  math- 
ematical processes  entailed  in  the  maximum  principle,  and  an  explicit 
formulation  for  the  Lagrangian  multiplier  functions.   This  explicit  for- 
mulation is  used  in  the  next  section  to  sho'w  the  necessary  conditions 
may  then  be  regarded  as  an  algorithm.   The  final  section  contains  a  brief 
summary  of  solution  methods. 

The  main  theoretical  development  of  the  dissertation  is  contained 
in  the  first  three  sections.   It  is  well  known  that  the  problems  encoun- 
tered in  the  calculus  of  variations  are  equivalent  to  the  optimization  of 

43 


44 


a  functional  (in  the  sense  of  mathematical  programming  problems)  under 
certain  restrictions  upon  the  variations.   A  good  exposition  of  this  is 
available  in  Luenberger  (1969).   With  this  equivalence  in  mind,  it  is 
noted  that  the  PMP  is  itself  worded  as  a  constrained  optimization  prob- 
lem.  I'Jhen  treated  with  what  is  normally  regarded  as  a  numerical  method, 
the  gradient  projection  method,  an  explicit  formulation  of  the  atten- 
dant Lagrangian  multipliers  is  obtained.   This  form  satisfies  all  of  the 
calculus  of  variation  necessary  conditions  and  allows  one  to  use  them 
in  a  most  straightforward  fashion.   As  a  result,  these  necessary  condi- 
tions may  be  directly  used  in  the  form  of  an  algorithm  to  obtain  a  solu- 
tion.  Furthermore,  it  is  believed  that  treating  the  PMP  as  a  mathemat- 
ical programming  problem  in  conjunction  with  the  gradient  projection 
method  helps  to  explain  the  effect  of  combined  control-state  constraints 
upon  the  maximum  principle. 

3.1   Problem  Statement  and 
Necessary  Conditions 

A  general  problem  which  represents  a  large  class  of  structural 

optimization  problems   is  treated  in  the  sequel.   The  functional 

J  =  /   L  (x,u)  dt  (3.1.1) 

0 

is  to  be  minimized  with  respect  to  the  control  _u(t)  where  the  state 
x(t)  must  satisfy  certain  boundary  conditions  and  a  differential  con- 
straint; in  addition,  an  inequality  constraint  involving  both  the  state 
and  control  must  be  satisfied.   For 
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u'^(t)  =  [u,  (t)  u„(t)  ...  u  (t)]  (3.1.2) 

—         1     z         m 

x'^'Ct)  =  [x^(t)    x^(t)    ...  x^(t)]  (3.1.3) 

the  subsidiary  conditions  to  minimizing  the  cost  function  J  are: 

X  =  l(x,u)  (3.1.4) 

Specified  Boundary  Conditions  on  x(t)  (3.1.5) 

<^^(-^,u)    <  0  Jl  =  l,...,q  (3.1.6) 

Terminal  time  t   is  considered  to  be  constant;  allowing  it  to  be 
r 

unspecified  requires  only  a  slight  modification  to  the  following 
derivation. 

This  problem  is  a  particular  form  of  a  very  general  one  treated 
by  Hestenes  (1966).   His  results  are  a  set  of  necessary  conditions  which 
must  be  satisfied  by  the  optimal  solution  and  include  the  maximum  prin- 
ciple.  To  obtain  the  necessary  conditions,  the  inequality  constraints 
are  converted  to  equality  constraints  in  the  manner  of  Valentine  (1937). 
These  constraints  and  the  differential  constraints  are  then  adjoined  to 
the  cost  function  via  Lagrangian  multiplier  functions  '^^^(t)  and  p.(t) 


respectively. 


<}>^(x,u)  +  s2(t)  =  0 


where  the  slack  variables   s  (t)  are  defined   such  that 


The  symbol  "="  denotes  "is  defined  by." 
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s^(t)  =  [-<f.j^(x,u)]^  >  0 


■J  =  /   L„(x,u)dC  +  /   p  (t)  [x-^(x,u)]  dt  + 
0  0 


/    U„(t)  [m.(x,u)  +  s2(t)]  dt 


Implied  summation  convention  is  used  whenever  a  vector  formulation  leads 
to  possible  ambiguities  in  later  developments.   Integrating  the  second 
integral  by  parts  gives  a  result  that  leads  to  the  variational  Hamiltonian. 


J  =  2.  21 


+  /  [^Q-l^l+   vij^*j^-2iP  +  ^^£s2]  dt 


Define: 


A  T 

H(x,u,£)  =  L  (x,u)  -  P  (t)  l(x,u) 


(3.1.7) 


the  variational  Hamiltonian,  and  H  which  will  include  terms  arising 
from  the  inequality  constraint. 


H   =  -  H  -  p^^^ 
or 

k'^  =  2.   (t)  .fU,u)  -  L  (x,u)  -  M At)    (}>  (x,u) 


Hence 


J   =  2.  21 


T. 


-  /    [h"  +2i  P  -  ^s2]  dt 
0     0 


VJith  the  exception  of  the  maximum  principle,  all  of  Hestenes'  necessary 
conditions  are  obtained  from  the  requirement  that  the  first  variation 
of  the  cost  function  vanish.   In  the  following,  "&k"   designates 
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"the  variation  of  2i"»  &  subscript  vector  designates  the  partial  deriva- 
tive with  respect  to  that  vector,  with  the  result  itself  a  column  vector. 
Thus , 


6J  =  p^  6x 


-  /    [6  J(K  -  i)  + 
0     0        - 


+  6uV  -  y,2s^6s^]  dt  =  0 


To  derive  the  PMP  requires  an  extensive  mathematical  development  and  is 
not  included  since  it  contributes  nothing  to  the  present  discussion. 
However,  the  necessary  conditions  are  listed  in  order  to  be  available 
for  later  reference. 


H  =  f(x,u) 
2_ 


£ 


=  -H 


X 


0  =  p.6x.      ,   t  =  0 
1   1 


0  =  p.6x. 
1   1 


■^ 


.   t  =  t. 


Specified  Boundary  Conditions  on  x.(t)  , 


0  =  H 


0  =  u^(t)  *^(x,u)  ,  yj(t)  >  0 


H(2inPT'  ii^PT'  -P)  ^  "(±nPT'  i^'  2) 


^PT'  -H3PT 


^PT' 


The  optimal  solution  must  satisfy  these  six  conditions  together  with  the 
Inequality  constraint  (3.1.6). 
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The  PMP  states  that  along  the  optimal  trajectory,  each  instant 
of  time  t,  state  JSr.pY^'^)  ^"^"^  adjoint  state  £(t) ,  treated  as  fixed,  the 
optimal  control  j-Up^^^^  ^^    that  admissible  control  which  minimizes  the 
variational  Hamiltonian.   In  the  present  context,  admissibility  requires 
that  _u(t)  be  piecewise  continuous,  the  set  of  admissible  controls  being 
denoted  by  fi.   Hence  the  PMP  indicates  that 

UQp^(t)  =  ARGMIN  [H(Xqp^,  ii  .  P ) ]  (3.1.8) 

ueQ, 

Notice  that  the  necessary  conditions  suggest  nothing  about  how 
a  solution  is  obtained,  but  merely  indicate  certain  functional  relation- 
ships that  must  be  satisfied.   However,  equation  (3.1.8)  seems  to  inti- 
mate that  solution  of  the  necessary  condition  of  PMP  involves  a  mathe- 
matical programming  problem. 

3.2   Mathematical  Programjning: 
Gradient  Projection  Method 

Having  shown  that  the  PMP  from  the  calculus  of  variations 
approach  to  an  optimization  problem  may  perhaps  be  related  to  a  mathe- 
matical programming  problem,  the  latter  will  be  discussed  in  general 
terms.   Consider  a  nonlinearly  constrained  optimization  problem 


x^p^  =  ARGMIN  [F(x)] 
X  e  fi 

subject  to 

g.(x)  <  0     j  =  1,  .  .  .  ,m 

where  P.   denotes  the  set  of  admissible  state  components  x.,  i  =  l,...,n, 
and  to  be  admissible  requires  only  the  satisfaction  of  the  m  inequalities. 
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Necessary  conditions  ^^?hich  2L-ipT  ^^^^   satisfy  are  given  in  the  Kuhn-Tucker 
theorem  : 

(i)   constraints  are  satisfied         §-(2inPT^  -  *-* 
(ii)   multipliers  exist  such  that  A.  >  0 

and  for  all   j  =  1,  .  .  .  ,in        '^i^i*'^PT^  "  ° 

m 

(iii)   and         V  F(x^„^)  +   E   X.  V   g . (x^„^)  =  0 

X  -OPT     .  -,  2      X  ^i  -OPT 
-  J=l  -'      -     -^ 

Observe  that  if  I   denotes  the  set  of  indices  associated  with  active 
constraints,  the  first  t^70  conditions  may  be  written  as 

j  e  I,   ^  g.(x)  =  0   and   A.  >  0 

j  «!  I,   ^  g.(x)  <  0   and   A.  =  0 


Fox  (1971,  pp.  168-176)  presents  a  very  readable  proof  of  this  theorem; 
a  more  mathematical  proof  using  vector  space  concepts  is  available  in 
Luenberger  (1969) . 

Many  methods  for  obtaining  a  numerical  solution  to  the  nonlinear 
programming  problem  described  by  the  first  two  equations  of  this  section 
have  been  developed.   The  gradient  projection  method  by  Rosen  (1960)  is 
used  frequently  in  structural  optimization.   Basic  to  the  method  is  the 
orthogonal  projection  of  the  cost  function  gradient  onto  a  subspace 
defined  by  the  normal  vectors  of  the  active  constraints.   An  inherent 
part  of  the  algorithm  is  the  concept  of  a  "feasible,"  "usable"  direction. 
Any  direction  d  is  feasible  if  an  increment  x   in  that  direction  improves 
the  cost  function,  i.e.,  decreases  F(x).  Direction  d  is  said  to  be 
usable  if  it  also  satisfies  the  constraints.   As  long  as  a  feasible. 


50 


usable  direction  exists,  the  cost  function  may  be  Improved.   A  constrained 
optimal  solution  2^p,T,  occurs  at  that  point  v/here  no  feasible  direction 
Is  also  usable,  i.e.,  any  attempt  to  improve  the  cost  violates  a  con- 
straint.  In  Appendix  B  these  concepts  are  used  in  a  concise  proof  of  the 
Kuhn-Tucker  conditions. 

Fox  (1971)  derives  the  matrix  P  which  projects  the  cost  function 
gradient  into  the  subspace  defined  by  vectors  normal  to  the  active  con- 
straints.  This  is  equivalent  to  subtracting  all  components  parallel  to 
vectors  that  are  normal  to  surfaces  of  active  constraints  from  the  nega- 
tive gradient  of  the  cost  function.   Recalling  the  definition  of  set  I  , 
consider  r  constraints  to  be  active  such  that 

I^  =  {a^,a2,...,a^} 

Define  a  vector  iv'hose  elements  are  the  corresponding  nonzero 
Lagrangian  multipliers,  and  another  vector  whose  elements  are  the  active 
constraint  functions 

a'^  =  [A    X         ...  A   ] 

1    2       r 

From  the  N^  vector,  a  matrix  N  is  introduced,  each  column  of  whi(;:h  is  the 
gradient  of  an  active  constraint.   Hence,  N  is  an  (n  >;  r)  matrix  where 


N  =  N   =  [N.  .] 


>i 


1  =  1 ,  .  .  .  ,  n 
and  >  (3.2.1) 


8N    9g 

N    =  — ^  =  -L 

ij    3x.    8x.  -^ 
i     1 


j  =  1, 
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With  these  definitions  of  _A  and  N,  the  third  Kuhn-Tucker  condition  can 
be  written  as 

V  F(x)  +  N  ^  =  0 

At  any  feasible  point  2£  where  g.(.x)    <   0,  the  direction  which  best 
improves  the  cost  function  is  the  negative  gradient  of  the  cost.   If 
those  directions  which  lead  to  constraint  violations  are  subtracted  from 
-V  F(x),  the  projection  matrix  P  is  obtained.   Directions  causing  a 

constraint  violation  are  specified  by  the  gradients  of  active  constraints, 

T 
i.e.,  the  columns  of  N  .   l-Hiat  is  required  of  S,  the  projection  of  the 

gradient,  is  that 

S  =  (-V  F(x))  -  N^  A  (3.2.2) 

where  A_  are  scalar  coefficients  to  be  determined  such  that  S  is  ortho- 

T 
gonal  to  each  column  of  N  ,  or 

T  T  -^ 
(N  )   S  =  0 
— 2S 

T 
Xfnen  the  matrix  equivalent  to  N   is  used  together  with  the  S  expres- 
sion (3.2.2),  the  result  is 

n'^(-V  F(x)  -  N  ^)  =  0 

such  that  the  ^  which  satisfies  this  orthogonality  condition  is: 

X   =  -    (N^'^K)"-'-  N^(V  F(x))  (3.2.3) 

—  2i    "' 

Unless  the  active  boundary  surface  normals  V  g . (x)  are  linearly  depen- 

2i  J 

T 
dent,  the  matrix  (N  N)  is  nonsingular.   Conversely,  if  this  matrix  is 
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singular  the  active  constraints  are  not  linearly  independent;  however, 
this  is  not  a  condition  encountered  in  most  real  cases. 

Substitution  of  the  _A  expression  into  the  S  equation  leads 
directly  to  the  projection  matrix  P: 

S  =  -P  V  F(x) 
X  — 


where 


P  =  I  -  N(n'^M)  "^  n"^  (3.2.4) 


where  I  is  the  identity  matrix.   The  direction  S  which  best  improves 
the  cost  is  given  in  terms  of  P,  where  P  and  N  are  given  by  (3.2.4) 
and  (3.2.1).   If  no  constraints  are  active  at  a  point  x>  then  N  is  a 
null  matrix,  P  reduces  to  an  identity  matrix,  and  the  direction  of  best 
improvement  is  coincident  V7ith  the  direction  of  steepest  descent. 

In  the  algorithm  associated  with  this  method  the  starting  point 
must  be  a  feasible  point  v;here  g .  (x)    <   0  for  all  j  =  l,...,m.   The 
design  then  proceeds  in  the  S  direction  until  the  solution  is  satisfied 
to  within  a  specified  position  tolerance  e.   Necessary  conditions 
generally  programmed  in  a  computer  program  are: 

|S.|^  e,         i=l,...,n 

A  .    >  0  .   j  e  I^ 
A  -    =  0  ,    j  ^  I^ 

It  is  readily  seen  that  for  S,  P,  and  \_   defined  as  above,  these  are 
completely  equivalent  to  the  Kuhn-Tucker  conditions. 
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3. 3   Gradient  Projection  Method  Applied 
to  the  Maximum  Principle 

Based  upon  the  preceding  discussion,  the  similarity  between  PMP 

and  the 'mathematical  programming  problem  can  be  discussed.   The  maximum 

principle  states  that  the  optimal  control  u^px  minimizes  the  variational 

Hamiltonian  with  respect  to  all  admissible  _u.   Or,  at  each  time 

0  <  t  <  t„,  ii^prp  minimizes  H(x,_u,p^)  with  respect  to  u  for  given  x   and  p^ 

and  where  u^prp  is  subject  to  constraints  ^q(^,u)    <   0,    H   =   l,...,q. 

Treating  this  as  a  mathematical  programming  problem,  the  following 

correspondences  are  recognized 

X  '^  JJ 

F(x)  ^   H(x,u,p) 

g.(x)  '^  <i>^(^y}^) 

A.  "^  vAt) 
3  ^ 

7  F(x)  -^  H 
X   —      JJ 

S  '\'  H         plus  constraints 
_u 

Contini.iing  to  identify  corresponding  quantities,  at  each  time  t,  let 

I.  denote  the  set  of  active  constraints,  taken  to  be  r  in  number. 
A 

^A  "    ^"l'    "">'    ■  ■  ■  '    "r^ 

Then  N"""   'v   <J)'^  =    [(J)        i>        ...(})      ] 

12  r 

a"^  ^^  y'^  =    [VJ      (t)    ti      (t)    ...    V      (t)] 
—         —  a.,  a.-,  a 

12  r 


K.i;=    [^..] 


L-      J    - 


j    =    1, . . . ,n 
j    =    1, . . . ,m 
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Furthermore,  define  (H  )   to  be  the  gradient  of  H  with  respect  to  u 
where  all  components  that  cause  a  constraint  violation  have  been  removed. 
Since  projection  matrix  P  removes  cost  function  gradient  components  that 

lead  to  constraint  violations,  consider  its  use  in  the  maximum  principle. 

T 
With  the  correspondent  to  N  identified  as  i>    ,  then 

~u 

P  =  I  -  1^(1  T*u^"^^  T  (3.3.1) 

—  u.  —     u 
and 

(H  )„  =  P  H  (3.3.2) 

u_  P      _u  ' 

From  the  Kuhn-Tucker  conditions,  this  implies  that  along  the 

optimal  trajectory  (t,>^p^,  iigp^) 

§  =  -  (H^)p  =0  (3.3.3) 

Similarly,  at  each  time 

p(t)  -    (±  ^  i^)"-^  1  T  "u  (3.3.4) 

ju   —     li   ~ 

from  which  it  follows 


u   —      u   — 


-^- 

VI   +  H 
-          u 

— 

0 

(H   + 

T 

= 

0 

(      - 

= 

0 

h" 

u 

= 

0 

(3.3.5) 


Or,  H  =  0  (3.3.6) 
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Hence  the  control  law  from  Hestenes'  necessary  conditions  can  be  derived 
from  the  PMP  condition  by  treating  it  as  a  nonlinearly  constrained  math- 
ematical programming  problem.   IJhile  using  the  gradient  projection  method 
in  the  derivation,  it  is  seen  that  equation  (3.3.5)  is  equivalent  to  the 
third  Kuhn-Tucker  condition.   The  second  Kuhn-Tucker  condition  is  iden- 
tical to  Hestenes'  necessary  condition  on  the  Lagrangian  multipliers  used 
to  adjoin  the  inequality  constraints  to  the  cost  function.   Satisfaction 
of  the  inequality  is  implied  by  requiring  the  first  Kuhn-Tucker  condition 
to  be  fulfilled,  where 

A.  >  0  -^       li^(t)    >   0  (3.3.7) 

Xjg.Cx)  =  0  ->yy(t)<J)^(x,u)  =  0  (3.3.8) 

g.(x)  <  0  ^     <('j^(x>u)  i  0  (3.3.9) 

Thus  by  treating  the  solution  of  the  necessary  conditions  of  the  max- 
imum principle  as  a  programming  problem  with  inequality  constraints, 
using  the  gradient  projection  matrix,  and  by  requiring  satisfaction  of 
the  Kuhn-Tucker  conditions,  an  explicit  formula  for  Hestenes' 
Lagrangian  multiplier  functions  has  been  derived.   It  is  further  demon- 
strated that  with  the  l-i„(t)  so  defined  satisfaction  of  the  extremum  con- 
trol law  condition  is  implied.   However,  before  this  treatment  can  be 
accepted  as  valid,  it  must  also  bs  sho^^m  that  the  system  of  canonical 
differential  equations  is  unchanged. 

Consider  the  state  system  equations 

X  =  H  =  f(x,u) 
-    2. 


56 


where 


h"  =  E'^(t)i(x,u)  -  Lq(x,u)  -  Ji'^CtHCx.u) 


It  is  obvious  that  the  explicit  form  of  p(t)  has  absolutely  no  effect 
upon  the  state  system  equation  expressed  in  canonical  form. 

Demonstrating  that  the  adjoint  system  equation  is  unchanged 
requires    the  method  described  by  Bryson  et  al.  (1964).   Consider  the 
general  problem  of  Section  3.1  again,  but  with  only  the  differential 
constraints  adjoined  to  the  cost  function,  i.e., 


T 
Min  {-J  =  p.  X. 


•^F     ^F 

+  /    (H  -  x'^£)dt}  (3.3.10) 

0     0 


subject  to:   4)n  (x»u)  ^0     £  =  l,...,q  (3.1.6) 

T 
where   H(x,u,£)  =  L  (x,u)  -  2.   (t)l(x,u)  (3.1.7) 

and     X  =  l(x,u)  (3.1.4) 

Again  let  I   denote  the  set  of  indices  associated  with  r  active 
A 

constraints  at  any  time  t 

(^^(x,u)  <  0      ,     ^  s^  ^A  ^  ^'ji^'"--'  "  ° 

and     ^     is  defined  as  before 
T 

The  problem  can  than  be  thought  of  as  minimizing  (3.3.10)  subject  to 
(t)(x,u)  =  0 
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I'Thile  on  the  constraint  surfaces  defined  by  this  equation  the  variations 
in  control  6u(t)  and  state  62S(^)  ^^^  ^^^   independent  but  instead  are 
related  through  the  subsidiary  requirement  that 


6^(x,u)  =  0 


or 


6x^  4)   (x J u^)  +  (5u  (})  (x jH^  =  0 


(3.3.11) 


This  imposes  a  restriction  to  the  admissible  variations.   For  cost  func- 
tion (3.3.10)  to  be  a  minimum,  it  is  necessary  that  its  first  variation 


vanish,  i.e. , 


6J  =  j2_  &x 


^F    *^F 
0     0 


T       T       T- 

6x  H  +  5u  H  -  6x  p 

—  X    —  u    —  -^ 


dt  =  0 


(3.3.12) 


It  has  already  been  shown  that 


A  T 

H   =  -  (H  +  u  (+.)   =  0 
u  u 


V7hich  will  be  used  to  advantage  shortly,  after  having  added  and  sub- 

T   T 
tracted  the  term  6u  (_U  ^)   from  the  integrand  of  (3.3.12). 


0  =  2.  i5  X 


0     0 


6x"'"(H^_  -  p)  +  Su'^H  + 


,  „  T,  T^.     ^  T.  T,. 
+  ou  (y  (h)       -    6u    {\!    q;) 
— ^u    — u 


dt 


Rearranging  terms  gives 


0  =  _£.  °2i 


^F     *^F 
0     0 


6x  (H  -  p)  + 
—   X   — 


+  6u'^(H  +  (u^i)    )  -  6u  (u  ±) 
—        u  ^u     — ^ 


dt 
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But, 


and 


Hence, 


(H  +  (y^f-)  =  (H  +  y"^^),  =  (-h")_,  =  0 
u     "...    —   —  -^u         u 


^F  'f 


0  =  ^  6x    +  / 

'o         0 


T       •       T  T 
6x  (H  -  2^)  -  6_u  4   y 


dt 


It  is  here  that  the  restrictions  imposed  by  the  active  constraints  are 
applied;  from  (3.3.11) 


„  T  ,T    .  T  ,T 
-Ou   6   =ox   di 
—  -^    —  —X 


such  that 


0  =  2^  6x 


0  =  £  6  X 


0     0 


^F    ^F 
+  f 


T       •       T   T 
62c  (H  _  -  p^)  +  62<   A_  y_ 

X  A 


dt 


T       •     T 

6x  (H  -  p  +  4   y) 

—   X   —    X  — 


dt 


0     0  L 

Applying  Euler's  lemma,  for  arbitrary  variations  in  the  state  v/hich 
satisfies  the  constraints. 


(\   -  P  +  ^  ii)  = 


=  0 


which  by  the  following  manipulations  is  shoxjn  to  be  the  adjoint  system 
equation  of  Hestenes. 


£  =  H^  +  (y  ^)^ 


(-H  -  y-'cf.) 


p  =  -  H 
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Thus,  the  explicit  formulation  for  ]j_(t)  obtained  by  applying  the  gradient 
projection  method  to  the  PMP  satisfies  all  the  necessary  conditions  of 
Hestenes. 

It  may  happen  that  in  some  cases  the  constraint  upon  control 
does  not  depend  upon  the  state.   It  can  be  shotim  that  the  p^(t)  explicit 
formulation  is  equally  valid  in  this  instance.   By  examination  of 
equations  (3.3.1)  through  (3.3.9)  it  can  be  verified  that  all  the  neces- 
sary conditions  except  the  adjoint  system  equation  are  satisfied.   To 
demonstrate  the  latter,  recall  that  when  on  a  constraint  boundary  the 
first  variation  of  both  the  cost  functional  and  the  constraint  function 
must  vanish.   That  is,  for 

^(u)  =  0 
both 

6J  =  0 
and 

6(|)  =  6u^  4)~  =  0  (3.3.13) 

To  derive  the  desired  equivalence,  the  same  term  must  be  added  and 
subtracted  from  the  integrand  of  &J   as  before,  again  arriving  at  the 
result 


T 

0  =  £  6x 


+  / 
0     0 


T       •       T   T 
&x    (H   -  P)  -  5u  ^ 


P 


dt 


IvTien  the  constraint  variation  (3.3.13)  is  introduced  into  this  last 
equation,  then  by  Euler's  lemma 

(H  -  p)  =  0 
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Since  it  was  stipulated  that  ^(u)  is  not  a  function  of  x,    the  equation 
may  be  written 

(H  +  \L^±)^   -  p  =  0      . 

(   -  H*  )   -  i  =  0 


£  =  H 

Thus,  the  expression  for  _p(t)  is  valid  when  the  constraint  inequality- 
depends  only  upon  the  control  ju(t)  . 

3.4  Maximum  Principle  Algorithm 

In  the  introduction  to  this  chapter  it  was  stated  that  the 
Lagrange  type  problem  from  the  calculus  of  variations  is  equivalent  to 
an  ordinary  mathematical  programming  problem  based  on  the  Kuhn-Tucker 
conditions.   Furthermore,  when  inequality  constraints  are  present  the 
necessary  conditions  are  equivalent  to  the  Kuhn-Tucker  conditions. 
It  was  demonstrated  in  the  preceding  section  that  if  the  PMP  is  itself 
treated  as  a  mathematical  programming  problem,  application  of  the 
gradient  projection  method  provides  an  explicit  solution  for  the 
Lagrangian  multipliers  associated  with  active  constraints.   This  explicit 
solution  for  y  (t)  also  satisfies  all  of  the  other  necessary  conditions 
for  an  optimal  solution.   The  ability  to  determine  y^Ct)  explicitly 
in  terms  of  parameters  and   functions  that  describe  the  problem  suggests 
the  possibility  of   converting   the  necessary  condtions  of  an  optimal 
solution  into  an  algorithm  for  obtaining  it. 
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Ensuing  discussion  of  the  algorithmic  form  of  the  necessary 

conditions  contains  the  implicit  assumption  that  all  equations  are  valid 

along  the  optimal  trajectory.   It  is  further  assumed  that  the  problem 

under  consideration  is  that  one  described  in  equations  (3.1.1  -  3.1.6). 

The  algorithm  requires  that  x(t)  and  £(t)  be  known  at  each  time 

0  <  t  <  t_  for  which  the  solution  procedure  is  as  follows. 
F 

(i)   Use  PMP  on  the  variational  Hamiltonian  to  determine  an  optimal 
control  _u  (t)  independent  of  constraints. 
u"(t)  =  ARGMIN  [H(x,u,2.)] 

Evaluating  the  inequality  constraints  with  u_  =  _u   reveals  which  of 
the  X.  =  l,...,q  constraints  are  active.   Let  r  denote  the  number 
of  active  constraints  and  I   the  set  of  indices  associated  with 
them. 

^A  "  ^"l'  ^2'  *  ■  •  '  "r-^ 

*j,(2£'U)  =0     ^  ^  ^A 

(f,^(x,u)  <  0     ^  ^  ^A 

From  this  the  vectors  whose  elements  are  the  nonzero  Lagrangian 

multipliers  and  corresponding  constraint  functions  are  defined, 

respectively,  at  the  instant  of  time  t. 

T 
y  (t)  =  [u    M    ...    n      ] 
-         a,   a       o. 
12       r 

12      r 
(ii)   Having  identified  which  of  the  q  constraints  are  active,  r 

components  of  lV,p^  are  specified  by  _(^(x,u)  =  0.   They  may  be 
solved  by  using  the  implicit  function  theorem,  which  requires 
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(j)   to  be  of  rank  r.   This  in  turn  requires  the  r  constraints 

— TJ 

which  are  active  at  point  x(t)  to  be  linearly  independent. 
To  determine  the  remaining  (m-r)  components  of  u^-jpr,,  requires 
that  \i_   be  known  at  time  t,  but 

ii(t)  =  -  (1  T  -^^'^  i  T  "u 

This  value  of  _m  is  used  to  determine  the  "constrained"  Hamil- 
tonian, 

h"  =  -  (H  +  )^±) 

A 

(iii)    With  the  nonzero  Lagrangian  multipliers  _p  known  and  H   conse- 
quently defined,  the  remaining  (m-r)  unknown  components  of 
u^   are  determined  from  the  control  law  for  the  constrained 
system,  i.e,,  from 

* 

H  =  0 

Once  u^p,p  is  completely  knovm,  the  adjoint  system  equations 
are  determined  by 

J, 

i  =  -  if 
2i 

The  process  outlined  above  then  allows  u^_,^  to  be  written  as 

— ur  i 

u^p^  =  ARGMIN  [H(x^p^,u,j,)] 

u  e  fi 

since  the  u   obtained  in  this  fashion  satisfies  if)  (xjjj)  ^  0  which  is  the 
only  requirement  for  being  admissible.   However  it  must  be  recalled 
that  these  equations  are  valid  along  the  optimal  trajectory;  it  remains 
to  be  shown  that  this  algorithm  may  be  employed  in  some  manner  to  obtain 
that  optimal  trajectory  and  to  demonstrate  their  satisfaction  along  it. 
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3.5   Solution  Methods 

Necessary  conditions  from  the  calculus  of  variations  provides 
a  Two  Point  Boundary  Value  Problem  (TPBVP)  to  be  solved,  which  is  in 
general,  nonlinear.   For  all  but  the  most  simple  problems  no  analytical 
solution  is  possible  and  if  any  solution  is  to  be  obtained  a  computer 
must  be  used  with  som.e  numerical  method.   A  discussion  of  the  available 
methods  and  their  relative  advantages /disadvantages   is  not  included  here 
due  to  the  availability  of  such  discussions  in  the  literature,  e.g., 
Bullock  (1966).   All  of  the  methods  involve  some  iterative  scheme,  and 
for  optimal  control  can  be  separated  into  t\J0  general  categories. 

(i)  Indirect  methods.   Schemes  which  require  an  intial  guess  of 
the  state's  solution:   In  these  methods  the  starting  point  is 
an  initial  guess  of  the  time  history  of  the  solution.   The  con- 
trol associated  with  the  solution  is  a  subsequent  calculation. 
Iteration  continues  until  the  state  satisfies  some  criterion 
connoting  convergence;  the  final  control  history  at  conver- 
gence is  the  optimal  control, 
(ii)  Direct  methods.   Schemes  which  require  an  initial  guess  of  the 
control  function:   The  starting  point  for  these  methods  is  an 
initial  guess  of  the  control  time  history.   For  this  class  of 
methods  the  state  associated  with  the  control  is  a  subsequent 
calculation.   Iteration  continues  until  the  control  satisfies 
some  convergence  criterion. 

The  method  of  quasilinearization  was  selected,  based  upon  the 
success  of  Boykin  and  Sierakowski  (1972)  in  applying  it  to  a  constrained 
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structural  optimization  problem.   Excellent  convergence  for  their 
problem,  the  capability  to  handle  nonlinear  systems,  and  the  avail- 
ability as  an  IBM  SHARE  program,  ABS  QUASI,  dictated  its  selection. 
In  the  application  to  the  examples  in  Chapters  IV  and  V  the  program 
required  no  modification.   As  a  result,  a  detailed  discussion  of  the 
method  of  quasilinearization  is  not  included. 

The  problem  discussed  in  preceding  sections  of  this  chapter 
falls  into  the  general  class  of  problems  that  QUASI  handles,  that  is, 

i  =  ^(Y,t) 
with  the  boundary  condition  of  the  form 
B^Y(O)  +  B^  Y(tp)  +  C^.  =  0 


where  t^,  square  matrices  B   and  B  ,  and  vector  C  ,  are  specified, 
constant  quantities.   The  specific  form  of  B  ,  B  and  C   depend  upon  the 
given  boundary  conditions.   As  described  in  algorithm  form 


In  terras  of  the  general  QUASI  nomenclature. 


Y  = 


R  J 


^(Y't)  = 


G^(x,£) 


Boundary  conditions  are  determined  by  those  specified  for  the  original 
system  and  by  the  necessary  conditions  outlined  in  the  first  section 
of  this  chapter. 

In  Chapter  VI  the  problem  treated  by  Boykin  and  Sierakowski  (1972) 
is  solved  by  the  gradient  projection  method  applied  to  a  finite  element 
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formulation  for  the  description  of  the  structural  system.   This  is 

a  method  of  the  second  kind  mentioned  above.   Results  of  the  two  methods 

are  compared. 


CHAPTER  IV 


CONSTRAINED  DESIGN  OF  A  CANTILEVER  BEAM 
BENDING  DUE  TO  ITS  OWN  WEIGHT 


4.0  Introduction 

A  structural  optimization  problem  has  been  selected  for  its 
simplicity  and  stated  as  an  optimal  control  problem.   The  maximum  prin- 
ciple is  applied,  giving  a  nonlinear  TPBVP  of  the  Mayer  type.   Among  the 
earliest  expository  papers  on  the  maximum  principle,  Rozono^r  (1959)  give; 
an  excellent  treatment  to  a  similar  type  of  problem;  his  technique  is 
used  to  obtain  both  the  variational  Hamiltonian  and  adjoint  variable 
boundary  conditions.   It  is  shoi^i  that  no  finite  solution  exists  for  the 
situation  of  unconstrained  control.   Numerical  solutions  for  constrained 
control  are  obtained  by  the  method  of  quasilinearization.   Constraints 
include  both  geometric  limitations  to  control  as  well  as  m.aximum  stress 
limits  that  become  mixed  constraints  depending  upon  both  state  and  con- 
trol variables. 

4 . 1  Problem  Statement 

A  cantilever  beam  of  variable  rectangular  cross  section  is  to  be 
designed  for  minimum  tip  deflection  due  solely  to  its  Ovvm  weight.   The 
m-itcrial  is  specified  to  the  extent  that  the  modulus  E  and  density  p  are 
constants.   Length  L  is  specified  but  the  design  variables,  height  h(x) 
and  ^^/idth  w(x),  may  be  chosen  independently  of  each  other,  subject  to 
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hard  constraints  upon  the  allowable  dimensions.   That  is 

a  <  w(x)  <  c 

(4.1.1) 
b  <  h(x)  <  d 

If  Y(x)  denotes  the  deflection  of  the  centerline,  the  problem  is: 
given  E,  p,  L,  and  the  constraints,  find  h (x)  and  w(x)  to  minimize 
Y(L).   The  particular  form  of  differential  constraints  to  be  satisfied 
will  be  derived  in  the  next  section. 

4.2   Structural  System 

Small  deflections  are  assumed  in  order  to  use  linear  Bernoulli- 
Euler  bending  theory.   Basic  conventions  assumed  for  this  example  are 
depicted  in  Figure  4.1;  with  these  conventions  the  governing  equation 
is  derived  using  standard  strength  of  materials  considerations.   The 
result  is 

EIg(x)Y"(x)  =  Mg(x)  (4.2.1) 

where 

L 
Mg(x)  =  y/   (t-x)w(t)  h(T)dT 

X 

Ij5(x)  =  Y^  w(x)  h3(x) 

and  Y  =  Pg-   Kinematic  boundary  conditions  to  be  satisfied  by  the  solu- 
tion of  (4.2.1)  are: 

Y(0)  =  0 
Y'  (0)  =  0 
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h(x) 


w(x) 


dW  =  Yw(x)h(x)dx 


/ 
/ 


/ 


/ 


Vg(x) 


M„  (x) 


Note:   Y(x)  is  centerline  deflection,  positive  downward. 


Figure  4.1   Structural  Conventions 
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Design  variables  and  the  related  constraints  (4.1.1)  are  put  into 
dimensionless  form  such  that 


^1  =  -^  -^       b/d  <  V,  <  1 
id  1 


w(x)        , 
v„  =  — —  ->   a/c  <  V-  <  1 
/     c  2 


and 


Ig(x)  =^cd3  v^v^ 

Replacing  the  independent  variable  with  a  dimensionless  equivalent,  and 
using  the  control  components  allows  the  governing  equation  to  be  put 
into  a  dimensionless  form.   For 

1   E   ,d,2  ^ 


3^-^(^)   uJu^Y  =  /   (x-t)  u^(T)u2(T)dT  (4.2.2) 

where 


yL2   "    -   -  t 


u.(t)  -  v.(x(t)) 

T'Jhen  constant  C   is  defined,  the  usual  kinematical  relationships  for 
a  beam  may  be  written  in  a  simple  dimensionless  form;  that  is,  let 

Cg  =  -y (y}  (Units  =  Length   ) 

X  =  C  Y  Deflection 

J-    ij 

x„  =  C  Y  Slope 

X  =  C  u^u„Y  Moment 


^4  =  S  ^  ^"1^2^^  ^^^^^ 


■,2 

u  u  =  C  (u|u  Y)  Load 

dt^ 
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These  state  component  definitions  are  used  with  the  natural  boundary 
conditions  to  obtain 

X3(l)  =0 

x^(l)  =  0 

From  (4.2.2),  the  state  component  definitions,  and  the  above  boundary- 
conditions  it  follows  that 

i^  =  x^  x^(0)  =  0 

^2  "  ""s^^I^a  ""2^^^  "  ° 

X  =  X,  X  (1)  =  0 

X,  =  u  u  ^4(1)  "=  0 

These  equations  and  boundary  conditions  are  used  in  the  following  section 
to  precisely  state  the  problem.   The  solution  and  results  are  given  later. 

4. 3   Unmodified  Application  of 
the  Maximum  Principle 

In  terms  of  the  state  variables  defined  in  the  preceding  section, 

the  problem  can  be  stated  with  more  mathematical  precision.   Find 

u^p^Ct)  =  ARGMIN  [x^(l)] 
_u  e  fi 

subject  to: 

(i)  differential  constraints  ii  =  _£_(.^,u) 

(ii)  kinem^itic  boundary  conditions  x  (0)  =  x  (0)  =  0 

natural  boundary  conditions  x„(l)  =  x,(l)  =  0 

(iii)  hard  geometric  constraints     b/d  <  u,(t)  <  1 

a/c  <  U2(t)  <  1 
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According  to  terminology  in  the  calculus  of  variations  this  is  a  Mayer 

type  problem.   Among  the  early  papers  concerning  the  PJ'IP,  RozonoSr  (1959) 

applies  the  PMP  to  a  similar  problem  giving  a  geometric  interpretation 

to  the  function  of  the  adjoint  variables. 

In  Rozonoar's  problem  the  cost  is  a  generalization  of  the 

ordinary  Mayer  problem,  in  the  sense  that  the  cost  function  is  a  linear 

combination  of  the  terminal  state  components.   It  can  be  shown  via  the 

calculus  of  variations  that  to  minimize 

T 
J  =  c^  x(tp) 

where  c   is  a  vector  of  prescribed  constants,  the  necessary  conditions  are 


— r 


H  =  /(t)f  (x,u,t) 


X  =  H   =  f 
-         2.       - 

T 
£  =  -H   =  -  f  i 


0  =  H 
u 


0  =  /(O)  x(0) 

0  =  /(tp)[c^  +2(tp)] 


H(^PX'  iioPT'  ^'^^  -  "*^^PT'  -'  ^'^-^ 


That  is. 


Min[J]  ->  Max  [H] 
u        u 
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If  we  consider  the  cantilever  beam  problem,  the  forms  in  the 
necessary  conditions  are 

c  =[10  0  0]'^ 


H  =  pj^x^  +  P2X3/U3U2  +  P3X^  +  P^Uj^u^ 


^F=l 


ry^r 


-\ 


f(x,u)  / 


X3/U3U2 


V"l"2   y 


from  which 


Pl  =  0 


Pj^Ci)  =  -1 


P2  =  -Pi 


P2(l)  =  0 


P3  =  ~P2^"l"2 


P3(0)  =  0 


P4  =  -P3 


p^(0)  =  0 


Adjoint  variables  p-,  (t)  and  p„(t)  can  be  integrated  by  inspect 


ion 


P^(t)  =  -1 


P3(t)  =  -  (1-t) 


0  <  t  i  1 


(4.3.1) 


such  that 


H  =  -  x^  -  (l-t)x3/u3u2  +  P3x^  +  P/^u^u^ 
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VJith  this  result,  the  necessary  condition  for  control  to  minimize  tip 
deflection  is 

H^^  =   3(l-t)x3/u;^U2  +  p^u^  =  0 

Hu2  =  (l-t)x3/uju2  +  p^u^  =  0 

At  first  this  appears  to  be  a  contradiction  since  the  two  equations  can 
be  satisfied  only  by  the  trivial  solution  because  they  have  equivalent 
forms, 

\  (4.3.2) 

P^uju2  =  -  (l-t)x3  J 

Further  examination,  however,  leads  to  the  conclusion  that  when  the 
control  is  completely  unconstrained  there  is  no  horizontal  tangent 
plane  to  the  surface  H  =  H(u  ,u„). 

\Jhen  the  geometric  constraints  to  the  control  are  included, 
a  constrained  minimum  may  exist.   If  such  is  the  case,  the  raaximum 
value  of  H  occurs  on  the  boundary  of  admissible  control  space.   To  that 
end,  PMP  is  em.ployed  along  the  control  space  boundary  to  determine 
u    (t)  at  each  time  t.   Before  detailing  this  procedure,  it  is  neces- 
sary to  first  consider  some  structural  aspects  of  the  problem. 

By  definition  the  control  components  are  positive,  which  in 
turn  implies 

Load:    Aa  >  ^  ^  ^4  "^  "l'^2 

Shear:    x   <  0  ^  ^a    ^   ^      '   x^CD  =  0 

Moment:   x^  >  0  ^  x^  <  0   ,   x^Cl)  =  0 
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Furthermore,  since  p„(t)  <  0   from  (4.3.1) 

P3  >  0  ^  P3  >  0  ,   p^CO)  =  0 

p^  <  0  ^  P4  <  0  ,   p^(0)  =  0 

This  exercise  makes  it  possible  to  use  the  information  of  the  sense  for 
X-  and  p,  to  simplify  the  search  for  u^p„  on  the  control  boundary.   By 
arranging  the  Hamiltonian  in  the  following  fashion 

H  =  -X2  +  p^x^  +  p^[u^U2  -  (l-t)x3/p^u3u2] 

it  is  observed  that  both  terms  in  the  bracketed  expression  are  positive. 
This  and  the  p.  outside  the  leading  bracket  allows  the  following 
equivalence: 

Max  [H]  ^   Min   [u-j^u^  -  (l-t)x3/p^u3u2] 
u  e  9U   u  e  8U 


or, 


where 


and 


u^    =  ARGMIN   [0(u)] 
u  c  3U 


■(u)  =  [u^u^  +  F2(t)/u3u^  ]  (4.3.3) 


F2(t)  =  -  (l-t)x3/p^  >  0 

At  each  position  t  along  the  beam,  the  state  and  adjoint  variables  must 
satisfy  the  appropriate  differential  equations,  and  _u  is  specified  by 
the  preceding  three  equations. 

Control  space  boundary  r)U  is  illustrated  in  Figure  4.2,  where 
the  u  axis  is  treated  as  the  ordinate  since  u  (t)  and  u  (t)  correspond 


b/d  <  u^  <  1 


a/c  <  u„  <  1 


75 


I-- 


U-. 


b/d-- 


a/c 


u^ 


Figure  A. 2  Admissible  Control  Space 
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to  the  height  and  width,  respectively,  of  the  cross  section  of  the 
beam  at  position  t. 

Along  the  constant  u^  edges  of  3 U,  let  u,  =  u  where  u  has  the 

-L  ice 

value  of  either  b/d  or  unity.   If  $ (u  ,u„)  has  a  minimum  point 

c  I  ' 

d-J    ^        ^      d2$ 

-, —  =  0     and     >  0 

^^•^9  A    2 


where 


'^"c'''2^  "  V2  ^   F2(t)/u3u2 

-^=u    -  F'(t)/u3u2 
du     c  c  2 


—  =        2F2(t)/u3u3 
du2  ^  2 


The  value  of  u„  which  satisfies  the  first  condition  is 

u„  =  F(t)/u2 
2         c 

Furthermore,  it  is  observed  that  only  one  extremum  of  $(u_)  exists  along 

u  =  u  and  that  it  is  a  minimum.   Hence,  either  ■i'Cu  ,u^)  has  a  minimum 
-L    c  c   2 

on  the  constant  u  edge  or  is  monotonlcally  decreasing/increasing.   If 

either 

*  k 

u„  <  a/c        or       1  <  u 

then  along  the  constant  u   edge,  H  has  its  maximum  value  at  a  corner  of 
the  rectangular  8U.   On  the  other  hand,  if 
a/c  <  u"  <  1 

then  H  has  its  maximum  value  on  the  line  u,  =  u  interior  to  the 

1    c 

endpoints. 
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Similarly,  along  the  u  edges  of  8U,  denote  u  =  u  where  u 
has  the  value  of  either  a/c  or  unity.   If  ^(u  ,u  )  has  a  minimum  point 

- —  =  0     and     — 2"  >  0 
du,  du 

1  1 


where 


$(u, ,u  )  =  u,u  +  F2(t)/u3u 
1   c     1  c  ±  c 


^^  u  -3F2(t)/u'^u 


du-  c  1  c 


du? 


12F2(t)/u5u 


It  is  observed  that  $(u, ,u  )  has  only  one  extremum  along  u„  =  u  , 

1   c  2    c 

it  is  a  minimum,  and  occurs  at  the  point  u^  =  u  where 

u*  =  +  {3F2(t)/u2}^ 
1  c 

Thus,  by  the  same  argument  posed  in  the  preceding  paragraph,  if  either 
u^  <  b/d     or     1  <  u 

then  along  the  constant  u   edge,  H  has  its  maximum  value  at  a  corner  of 
the  rectangular  3U.   VJlierever 

b/d  <  u*  <  1 

H  has  its  maximum  value  on  the  line  u„  =  u   interior  to  the  endpoints. 

2     c 
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On  the  basis  of  these  arguments,  the  following  system  was 
solved  by  the  method  of  quasilinearization: 
i^  =  x^  x^(0)  =  0 

x-  =  x„/u|u2        x„(0)  =  0 

x„  =  X,  X  (1)  =  0 

3    4  3 

P3  =  -P2^nlu^     ■       ■'   P3(0)  =  0 
P4  =  -P3  P^CO)  =  0 

-"-OPT  "  ARGMIN  [u^u^  -  (l-t)x3/p^u3u2] 

The  beam  is  represented  by  100  intervals  composing  the  range  0  <  t  <  1, 
which  is  separated  by  101  "mesh  points."  An  initial  guess  of  the 
solution  x(t)  and  p(t)  is  chosen;  it  is  selected  to  satisfy  the  bound- 
ary conditions.   This  guess  is  not  a  solution  and  does  not  satisfy  the 
differential  equations.   The  x  and  2_   equations  are  linearized  about  the 
Initial  guess,  then  the  resulting  linear  TPBVP  is  solved  to  obtain  new 
x(t)  and  p(t)  functions  which  more  closely  satisfy  the  differential 
equations.   At  each  time  t  corresponding  to  a  mesh  point,  H  is  numer- 
ically evaluated  along  each  of  the  four  straight  line  segments  com- 
posing 3U  to  determine  Uv^pr,,-   The  point  (u,  ,u  )  on  3U  which  gives 
H(ji;x,2,t)  its  maximum  value  is  u,,.prp-   This  process  is  repeated  until 
the  x(t)  and  £(t)  iterate  satisfies  the  differential  equations  to 
within  a  specified  tolerance.   The  equations  necessary  to  use  the  IBM 
program  available  are  given  in  Appendix  C,  in  the  form  of  a  subroutine 
listing. 
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4. 4   Results:   Geometric  Control  Constraints 

For  the  most  part,  no  major  difficulties  were  encountered  in 
using  quasilinearization  to  obtain  a  solution  to  the  sixth  order  sys- 
tem derived  in  the  previous  section.   Certain  parameter  values  did 
engender  numerical  instability.   These  cases,  the  source  of  the  diffi- 
culty, and  its  circumvention  are  discussed  in  Chapter  VII.   Moreover, 
all  calculations  were  done  in  double  precision  as  necessitated  by  matrix 
inversion  accuracy  requirements. 

Parameter  values  selected  to  illustrate  the  solution  method  are: 


u   :   b/d  =  0.25 


u   :    a/c  =  0.20 


(4.4.1) 


The  measure  of  error  of  satisfaction  of  the  differential  equations  in 
the  TPBVP  is  in  terms  of  the  general  system 

dY. 

-—  =f.(Y,x)     Y  =  Y.(x),      i  =  l,...,n 
dx     1  —  —1 

ERROR  -  Max  I dY .  -  f.(Y,x)dx|  (4.4.2) 

1 

Deflection  of  a  uniform  beam  due  to  its  own  weight  was  used  to  infer 

an  initial  guess  which  satisfies  all  boundary  conditions: 

x^(t)  =   t^^ 


0  <  t  <  1  (4.4.3) 


X2(t) 

= 

t3 

X3(t) 

= 

1    - 

t 

x^(t) 

= 

-1   + 

t 

P3(t)  =   t2 
P^(t)  =  -t3 
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With  these  specified  parameter  values  and  initial  guess  of  the  solu- 
tion, the  program  converged  to  a  solution  in  five  iterations.   From 
this  run  a  tolerance  was  selected  for  all  subsequent  cases;  the  follow- 
ing tabulation  provides  the  data  used  in  its  selection: 


teration 

ERROR 

Tip  Deflection  (Cos 

1 

.2028 

.7387749327 

X  10"-^ 

2 

.1704 

.3192152426 

X  lO"^ 

3 

.6533 

X 

lO" 

-1 

.2847993812 

X  10~^ 

4 

.3031 

X 

lO" 

-5 

.2853731846 

X  10~^ 

5 

.1129 

X 

lO" 

■10 

.2853719983 

X  10~^ 

It  is  seen  from  these  tabular  data  that  there  is  little  improvement  in 
cost  (tip  deflection)  as  a  result  of  the  fifth  iteration.   For  this 
reason  a  value  for  the  tolerance  was  selected  as  0.5  x  10   which  corre- 
sponds to  about  six  significant  digits  in  the  cost  functional. 

Recall  from  the  previous  section  that  no  unconstrained  minimum 
exists.   With  the  control  bounds  included,  the  intuitive  solution  is 
one  in  v/hich  the  cross-sectional  area  is  maximum  near  the  root,  and 
reduces  to  a  minimum  at  the  tip.   Recalling  that  for  the  optimal  control, 

n    p^==Max      [H]     =     Min      [0] 
_u  e  ii  11  e  i2 

A  sequence  of  illustrations  in  Figure  4.3  demonstrates  the  location  of 
u    on  oU  for  several  stations  along  the  beam.   Constant  contours  of 
^(u^)  are  plotted  on  the  admissible  control  space  at  five  distinct  posi- 
tions.  If  an  extremal  point  exists  interior  to  3U  some  lines  of  constant 
*(u^)  contours  must  be  closed  curves  in  jj-space.   This  is  impossible  for 
this  example. 
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x/L  -  .1 


u. 


0  1 

"2 


x/L  =  .3 


0   1 


x/L  =  .5 


0  1 


x/L  =  .7 


x/L    =    1. 

n 

i\  w 

\^ 

(L^ — _^;:^ 

- 

O  Minimum  $(_u)  ^-  -^0?! 


Lj  Maximum  §(u) 


Direction  of  Increasing  ^Cu) 


Figure  4.3   Contour  Plots  of  i>iu)    at  Various  Stations 
Along  the  Beam 
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The  first  illustration  is  for  the  t  =  0.1  cross  section,  near 
the  root  of  the  beam.   Since  u^p™  occurs  at  the  point  of  minimum  'J'Cii), 
the  optimum  value  for  both  u,  and  u„  is  unity,   the  maximum  allowable 
dimensions  for  both  height  and  width.   Constant  contour  lines  indicate 
that  ^  (jj)  is  mathematically  decreasing  in  either  direction  of  ^-space. 
Lines  of  constant  '^(u)    are  also  plotted  for  the  cross  section  of 
t  =  0.3.   The  optimum  control  has  the  maximum  admissible  value  for 
height  u  but  u„  has  a  value  somewhat  less  than  unity.   However,  there 
are  still  no  contour  lines  which  are  closed  curves. 

At  the  midpoint  cross  section  the  minimum  $(_u)  point  occurs 
at  u  =  1  and  u  =  a/c.   Although  the  surface  <5(u_)  forms  a  scoop-like 
shape,  there  are  still  no  closed  curve  contours,  and  hence  no  extremal 
interior  to  admissible  u.-space.   The  next  cross  section  at  which  ^(_u) 
is  displayed  occurs  at  t  =  0.7.   On  this  section,  u„  is  still  at  its 
lower  bound  but  u   is  no  longer  at  the  maximum  allowable  value  of  unity 
as  shown  in  Figure  4.3.   In  the  last  of  the  sequence,  '!'(u.)  contours  for 
the  cross  section  at  the  tip  of  the  beam  are  displayed.   The  point  of 
m^inimum  '^(u)    occurs  v/here  both  components  of  control  have  their  minimum 
allowable  values.   Again,  no  contour  lines  of  constant  v(jj)  form  a 
closed  curve  indicating  the  existence  of  an  interior  extremal  point. 

This  sequence  of  illustrations  indicates  two  things.   First, 
the  lack  of  closed  curve  contour  lines  of  ^(u)    verifies  that  u^p™  exists 
on  9U.   With  further  study  it  may  be  possible  to  obtain  some  condition 
on  H   =0  which  implies  the  equations  corresponding  to  (4.3.2)  can 
never  yield  a  finite,  unconstrained  optimum.   Such  a  condition  would 
define  the  class  of  structures  whose  unconstrained  solution  is  the 
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"zero  volume  solution"  frequently  described  in  the  literature  on  struc- 
tural optimization.   Secondly,  at  t  =  0  the  point  u^p„  occurs  at 

T 
u  =  [1,1]  ,  the  point  of  maximum  cross-sectional  area;  as  t  increases 

from  zero  to  one,  the  point  iVvprr.  moves  along  the  u  =  1  boundary  of  3U 

to  the  u  lower  bound,  and  then  down  the  u  =  a/c  boundary  of  9U  to  u 

lower  bound.   By  the  time  t  =  1  the  optimal  cross-sectional  area  is  the 

minimum  allowable  area.   As  a  result  of  the  prescribed  form  of  9U, 

if  _n   (t)  follows  this  particular  path  as  t  increases  from  zero  to  one, 

each  component  of  u,,p„(t)  has  its  o^^m  distinct  region  of  transition. 

That  is,  at  any  value  of  t,  if  b/d  <  u,.  <  1,  then  u  must  be  on  either 

its  upper  or  lower  bound.    Conversely,  if  u^  is  in  transition  where 

a/c  <  u„  <  1,  then  u^  must  be  on  one  of  its  bounds. 

This  effect  is  seen  most  clearly  in  Figure  4.4,  where  u,,p,^  is 

displayed  for  the  example  case  parameter  values  specified  by  (4.4.1). 

The  profiles  are  displayed  on  a  two-view  drawing  as  a  plan-form  of  the 

beam  might  appear.   State  components  corresponding  to  this  beam  are 

shown  in  Figure  4.5,  representing  dimensionless  deflection,  slope, 

moment,  and  shear,  respectively.   As  observed  in  Figure  4.4,  there  are 

five  distinct  regions  of  the  beam: 

(i)   0  <  t  <  .25 


^1=  ^ 
U2=l 

(ii)  .25  <  t  <  .52 


^=^ 


controls  onupper  bounds 


a/c  <  u   <  1 


u-  transition 


%U2 
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1.0 


x/L 


a/c  =  .20 


TOP  VIEI^? 


SIDE  VIEW 


.5 


b/d  =  .25 


-gu^ 


x/L 


Figure  4.4   Plan-form  Views  of  Optimal  Solution  for 
b/d  =  .20  and   a/c  =  .25 
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Figure  4.5   State  Coraponents  of  Optimal  Solution  for 
b/d  =  .20  and  a/c  =  .25 
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(iii)   .52  <  t  <  .59 


"1=1 


u  =  a/c 


(iv)   .59  <  t  <  .90 


b/d  <  u^  <  1 

u  =  a/c 


(v)   .90  <  t  <  1.00 


u^  =  b/d 
u^  =  a/c 


on  upper  bound 
on  lower  bound 


u^  transition 


controls  on  lower  bounds 


The  curves  that  show  the  intercept  locations  as  a  function  of 
parameter  values  b/d  have  been  called  "correlation  curves"  in  earlier 
studies.   Ifhen  the  xv^idth  is  allowed  to  vary  also,  the  second  parameter 
a/c  is  introduced.   For  the  sake  of  comparison  to  previous  studies,  the 
intercept/correlation  curves  are  plotted  as  dependent  upon  b/d  and 
parametric  in  a/c.   However,  it  would  be  just  as  correct  to  do  the 
opposite. 

Intercept  location  curves  described  above  are  sho'.ra  in 
Figure  A. 6.   The  heav^'  black  curve  is  the  case  where  a/c  =  1,  a  beam 
of  constant  uniform  width — the  case  cited  from  earlier  literature. 
Another  special  case  is  represented  by  dashed  lines,  corresponding  to 
a/c  =  0  which  is  the  case  corresponding  to  a  minimum  allowable  thick- 
ness equal  to  zero.   Dashed  lines  are  used  because  these  data  are  an 
extrapolation:   convergence  problems  encountered  for  parameter  values 
less  than  0.1  prevented  obtaining  numerical  results. 
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A  discussion  on  the  convergence  difficulties  experienced  by 
the  quasilinearization  algorithm  for  parameter  values  approaching  zero 
is  presented  in  the  chapter  on  numerical  instabilities.   In  that  dis- 
cussion, isolation  of  the  source  of  difficulty  is  reported;  it  is  pos- 
sible that  this  difficulty  may  be  a  general  result  applicable  to  all 
problems  to  be  solved  by  the  method  of  quasilinearization.   A  solu- 
tion for  this  case  is  later  obtained  by  finite  element  techniques. 
Note  that  since  0  <  a/c  <  1  these  two  cases  represent  limits  to  the 
solutions  of  the  problem.   In  addition,  if  the  four  intercept  locations 

are  plotted  versus  a/c  and  parametric  in  b/d,  curves  r   and  r   appear 

a      c 

as  "horizontal  vees"  with  r.  and  r,  lines  that  are  nearly  parallel. 

b      d 

It  is  interesting  to  note  from  the  figure  depicting  the  solu- 
tion of  this  case  as  a  plan-form,  that  in  the  central  region  of  the 
beam,  the  height  is  greater  than  the  width.   This  result  can  be  antic- 
ipated since  such  a  configuration  gives  a  greater  bending  resistance 
per  unit  weight. 

With  further  reference  to  Figure  4.4,  the  transition  of  u^ (t) 
is  seen  to  be  almost  a  linear  taper,  whereas  the  u-(t)  transition 
exhibits  a  much  more  pronounced  curvature.   To  generalize  from  this 
specific  case  of  given  values  of  b/d  and  a/c  to  arbitrary  values 
requires  the  introduction  of  four  quantities  characterizing  the  solu- 
tion.  These  quantities  are  the  values  of  t  at  the  points  where  the 
transitions  intercept  the  bounds  on  u  and  u  ;  since  t  represents  a 
normalized  position  x/L,  these  quantities  can  be  thought  of  as  an 
intercept  location  expressed  as  percent  of  the  beam's  length.   They 
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are  defined  with  reference  to  the  five  distinct  regions  of  the  beam 
previously  given,  where 

r  designates  u, (t)  intercept  with  lower  bound 
r  designates  u^ (t)  intercept  with  upper  bound 
r   designates  u„(t)  intercept  with  lower  bound 

3.  Z. 

r   designates  u  (t)  intercept  with  upper  bound 

such  that  the  five  regions  are: 

(i)   0  <  t  <  r      control  on  upper  bounds 
-   -   c 

(ii)  r   <  t  <  r      u„  transition 
c        a      2 

(iii)  r   <  t  <  r,     control  on  upper/lower  bounds 
a  -   -   d 

(iv)  r^  <  t  <  r,      Ut  transition 
d        b      i 

(v)  r.  <  t  <  1.0    control  on  lower  bounds 
b  -   - 

4.5   Inequality  Stress  Constraints 

This  section  treats  an  inequality  limit  to  allowable  normal 
and  shear  stresses  associated  with  bending.   Using  the  ordinary 
strength  of  materials  formulations  it  can  be  shown  that  for  the  rec- 
tangular cross  section  these  constraints  take  the  form 

h(x)  <  a. 


and 


2  I^(x)  ^   '    -   MAX 


1  ^B^^'^ 

8  ITW  ^'^"^  -  ^MAX 


When   diraensionless   quantities   are   introduced 


6  Y  ^3/^1^2  -"  SlAX 
-  I  YL   \/"i^2  -<  "max 
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the  inequalities  may  be  vn-itten  in  the  required  form  for  mixed  con- 
straints, i.e.,  as  a  function  of  both  control  and  state  components: 

(f)^(x,u)  =  x^/uju^  -  ^0  -  °  (4.5.1) 

^2^—'—^^    ^  ~   ^4^"l^2  ~  ^0  -  ^  (4.5.2) 


where 


1  ^MAX  ,d, 
°0  =  6  "^  ^"f^ 


2  '^MAX 


0   3  yL 

The  tv7o  stress  constraints  place  restrictions  upon  the  minimum 
cross-sectional  dimensions  to  keep  the  normal  and  shear  stresses  less 
than  prescribed  values.   Specifically,  from  the  constraints  (4.5.1) 
and  (4.5.2),  two  control  inequalities  must  be  satisfied  at  each  station 
t,  and  these  inequalities  depend  upon  the  state  of  the  structural  sys- 
tem.  The  inequalities  are: 

''l"2  -  "  ''4^^0   '      ^4^^)  <  0 
from  which  can  be  derived  boundary  arcs  in  _u-space: 


(4.5.3) 


"lT^^^2^  =  -  ^4/^0^2  -  ° 


Both  of  these  boundary  arcs  are  hyperbolas  restricted  to  the  first 
quadrant  of  _u-space.   Depending  upon  the  location  of  the  arcs,  vis-a-vis 
the  rectangular  9U,  inclusion  of  stress  constraints  has  one  of  three 
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effects  in  determining  what  _u  is  admissible.   At  any  station  t  in  some 
structural  state  x(t) , 

(i)   if  a    ,T      is  too  small  the  stress  boundary  arc  lies  entirely 
above  rectangular  8U;  all  geometrically  admissible  u   violate 
the  stress  constraints, 
(ii)   if  cr  ,T   is  too  large  the  stress  boundary  arc  lies  entirely 
below  rectangular  9U;  all  geometrically  admissible  u   satisfy 
the  stress  constraints. 
(iii)   for  some  range  of  o„,t   the  stress  boundary  arc  divides  the 

rectangular  3U  into  two  regions:  the  upper  region  consists  of 
geometrically  admissible  u  that  satisfy  the  stress  constraint, 
u  in  the  lower  region  are  geometrically  admissible  but  violate 
the  stress  constraint. 

Inclusion  of  stress  constraints  alters  the  admissible  control 

space  from  the  rectangular  shape  previously  considered  to  a  shape  that 

may  contain  a  stress  boundary  arc  as  part  of  its  boundary.   Consider  the 

normal  stress  boundary  arc  specified  by  (4.5.3)  to  be  a  part  of  (fU. 

Then  to  find  u^p„  in  the  manner  outlined  in  Section  4.3,  $(u)  must  be 

evaluated  alonp  u,   =  u,  (u„).   If  a  minimum  exists  along  the  orthogonal 
la     la  Z 

projection  of  u   (u„)  on  the  ^(u)    surface,  then 

^80(u) 


Min 


9(u) 


la 


3u^ 
32$(u) 


=  0 


>  0 


V.  8u2 
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Along  the  normal  stress  boundary  arc  u   (u  ) , 
u2u2  =  X3/aQ  ->   u^  =  X3/aQu2 

Substituting  for  u   in'<J'(_u)  it  follows  from  (4.3.3)  that 

$(u)l    =  (xJa     +   F2(t)a  )/u 
—  u,      3   0         U    1 

Since  the  expression  in  parentheses  is  positive,  it  is  readily  observed 
that  the  value  u  =  +  »  minimizes  ^(u)    along  u   (u  ) .   From  (4.5.3),  the 
value  of  u„  at  this  point  in  _u-space  is  zero. 

The  same  result  exists  for  <I>(_u)  evaluated  along  u   ;  on  the 
shear  stress  boundary  arc 

Substituting  for  u  in   (u.)  it  follows  from  (4.3.3)  that 


^(u) 


-1 


-,9 


u   =  ^-^^'^O^  -'   ^-^4/'^"0^    F^(t)/u^ 

It 


Since  both  the  expression  in  parentheses  and  F^(t)  are  positive,  it  can 
be  argued  as  above  that  the  point  of  minimum  ^(_u)  along  u   -:(u   )    occurs 
at  the  point  (u,  _=  +  ">,  u  =  0)  . 

1  L  / 

At  any  cross  section  there  exists  a  single  minimum  of  0(u)  along 
either  stress  constraint  boundary  arc.   Since  these  minima  occur  at  the 
point  u:  (+",0),  as  one  proceeds  along  either  boundary  in  the  direction 
of  increasing  u.  ,  <?(ja)  is  a  monotonically  decreasing  function.   For  this 
reason,  along  the  stress  constraint  boundary  arc,  it  is  necessary  to 
e'/aluate  ^(u^)  at  only  one  point.   That  point  is  the  intersection  of  the 
stress  boundary  constraint  arc  with  the  rectangular  boundary  having  the 
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larger  value  of  u, .   It  is  possible  that  the  two  stress  boundaries 
intersect,  with  the  coordinates  of  that  point  given  by 


r^ 

ro] 

"it  = 

\-" 

,~x,  j 

\   0^ 

2 

(%] 

rv 

u        - 

21 

.  y-^ . 

I      T_  , 

V  3' 

^      0^ 

Even  though  this  adds  a  complication  to  the  boundary  of  what  is  admis- 
sible u,    the  monotonically  decreasing  property  of  both  stress  boundary 
arcs  still  results  in  the  necessity  of  evaluating  $(_u)  at  only  a  single 
point  along  the  stress  boundary  arcs.   To  illustrate  this.  Figure  4.7 
depicts  an  admissible  control  region  determined  by  a  combination  of 
geometric  constraints  and  a  composite  stress  constraint  boundary, 
together  with  the  point  at  which  <J>(u)  is  a  minimum  along  the  composite 
stress  boundary.   With  continued  reference  to  the  figure,  u^prj,  is  that 
point  along  line  ABC  which  minimizes  <i>(u). 

4 . 6   Results:   Stress  Constraints  Included 

For  the  purpose  of  comparison,  the  parameter  values  used  in  the 
case  with  only  geometric  constraints  were  also  used  in  the  cases  where 
stress  inequality  constraints  are  present.   The  results  are  quite  simply 
stated:   neither  stress  constraint  is  ever  active.   In  general  it  was 
found  that  the  critical  section  occurs  at  the  root.   At  that  section 
the  greatest  portion  of  geometrically  admissible  control  space  is  pro- 
hibited by  the  stress  constraints.   The  further  along  the  beam  towards 
the  tip,  the  less  the  prohibited  geometrically  admissible  control  space. 
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Both  x^  and  -x,  decrease  monotonically  to  zero  at  t  =  0,  and  is  the 
3       4 

reason  for  this  phenomenon. 

Each  of  the  parameters  a„  and  t_  thus  has  three  important 
critical  values.   For  the  first  pair  of  values  the  maximum  allowable 
stresses  are  sufficiently  large  that  all  geometrically  admissible 
control  contained  by  the  rectangular  9U  satisfy  the  stress  constraints. 
In  this  case 

a  =  9.128 

T  =  8.658 

and  the  stress  constraints  intersect  the  lower  left-hand  corner  of  8U 
at  t  =  0  and  lie  below  9U  for  all  t  >  0.   For  the  range  of  values 

.5705  t  ^0  1   9.128 
2.1645  <  Tp  <  8.658 

the  stress  boundary  arcs  intercept  the  u^  =  • 2  portion  of  rectangular 
3U.   If  they  are  exactly  equal  to  the  lower  values,  respectively,  both 
arcs  pass  through  the  point  (1.,.2)  at  t  =  0,  and  lie  below  that  point 
for  all  t  >  0.   However,  these  values  are  not  as  important  as  either 
the  first  pair  or  the  next.   For  the  range  of  values 

.1141  <  Oq  <   -5705 

.4329  <  Tq  <  2.1645 

both  arcs  intercept  the  u  =  1  portion  of  rectangular  au.   When  equal  to 
the  lower  values  both  arcs  pass  through  the  upper  right-hand  corner  of 
3U  at  t  =  0.   If  the  values  of 
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a  <  ,1141 
T  <  .4329 

then  at  t  =  0  there  is  no  control  which  satisfies  both  geometric 
constraints  and  stress  constraints.   The  problem  has  inconsistent 
constraints  such  that  no  solution  is  possible. 

A  sequence  of  plots  showing  the  admissible  control  space  is 
shown  in  Figure  4.8  for  six  separate  stations  along  the  beam.   The 
value  of  o   is  .1141  for  this  case  such  that  in  section  t  =  0  the  set 
of  admissible  controls  is  the  single  point  u:    (1,  1).   That  portion  of 
the  geometrically  admissible  control  space  disallowed  by  the  stress 
constraint  is  indicated  by  dashed  lines.    Additionally,  the  point  cor- 
responding to  u^p,p  is  included  for  each  section.   By  observing  the 
sequence  of  plots  as  t  increases,  both  the  stress  boundary  and  u^prp  are 
seen  to  move  closer  to  the  origin  of  u.-space.   However,  u^prp  never  over- 
takes the  stress  boundary,  thus  the   normal  stress  boundary  never 
becomes  active.   Figure  4.9  depicts  the  same  case  for  the  shearing 
stress  constraint. 

Two  additional  figures  are  included  v;ith  the  intermediate 
critical  values  of  o      and  x  .   The  normal  stress  constraint  case 
a      --   0.5705  is  shoxm  in  Figure  4.10;  shear  stress  constraint  case 
T   =  2.164b,  in  Figure  4.11.   No  illustrations  are  provided  for  the 
cases  associated  with  the  largest  critical  values  since  for  t  >  0  all 
geometrically  admissible  controls  satisfy  the  stress  constraints. 
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CHAPTER  V 


CONSTRAINED  DESIGN  FOR  AN  OPTIMAL 
EIGENVALUE  PROBLEM 


5. 0  Introduction 

A  column  is  to  be  designed  for  maximum  buckling  load  in  order 
to  demonstrate  (i)  an  application  of  the  theory  developed  in  Chapter  III, 
and  (ii)  an  optimal  eigenvalue  problem.   Difficulties  pertinent  to  this 
class  of  problems  are  discussed  as  well  as  a  method  to  obtain  a  solu- 
tion by  using  a  theorem  concerning  self-adjoint  problems.   Use  of  this 
theorem  with  the  maximum  principle  overcomes  the  difficulty  associated 
with  the  simultaneous  dual  optimization  required  for  problems  where  the 
cost  functional  is   an  eigenvalue  expressed  as  a  Rayleigh  quotient. 
Numerical  solutions  to  the  combined  set  of  state  and  adjoint  variables 
are  obtained  by  the  method  of  quasilinearization.   Characteristics  of 
the  solutions  for  a  variety  of  control  bounds  are  presented,  along  with 
a  discussion  of  attempts  to  also  include  a  mixed,  inequality  constraint 
related  to  stress. 

5. 1  Problem  Statement 

A  vertical  column,  fixed  at  the  base  and  free  at  the  tip,  is  to 
be  designed  such  that  the  buckling  load  is  maximized.   Length  L  is 
a  specified  constant  as  is  the  total  weight  W.   The  dependence  of  cross- 
sectional  area  A(x),  material  modulus  E(x),  and  density  P (x) ,  upon 
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position  X  along  the  colunm  is  to  be  such  that  the  vertical  end  load  P 
at  buckling  is  maximized.   Weight  is  not  to  be  neglected  and  all  cross- 
sectional  shapes  are  similar. 

In  more  specific  terms,  a  given  weight  of  material  is  to  be 
arranged  with  variable  properties  and  geometry  such  that  the  result- 
ing mass  distribution  M(x)  and  stiffness  distribution  S(x)  maximizes 
the  vertical  tip  load  the  column  can  support  when  the  effect  of  weight 
is  included.   Upper  and  lower  bounds  to  all  design  variables 

A^  <  A(x)  <  Ay 

E   <  E(x)  <  Ey  (5.1.1) 

\  1   P(>^)  <  ^U 

are  also  specified.   Governing  equations  for  the  structural  system  are 
derived  in  the  next  section. 

5.2  Structural  System 

In  the  subsequent  analysis,  small  deflection  bending  theory  is 
assumed:   the  result  is  a  linear  differential  equation  V7ith  variable 
coefficients.   Without  this  assumption,  the  governing  equation  is  that 
of  the  nonlinear  elastica  problem.   Additionally,  axial  compression 
effects  are  neglected  as  second  order.   Let  Y(x)  denote  the  deflection 
of  the  centerline  at  position  x  along  the  beam  as  depicted  in  Figure  5.1 
Since  the  system  is  conservative,  the  governing  equations  can  be 
derived  quite  simply  by  energy  techniques.   In  terms  of  stiffness  S (x) 
and  mass  per  unit  length,  the  equation  and  boundary  conditions  are: 
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V^(x) 


MgCx) 


Figure  5.1   Structural  Conventions 
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L 
(S(x)Y")   +  PY"  +  g(Y'  /    M(C)d5)   =  0  (5.2.1) 

C=x 

Y(x)  =  Y'  (x)  =  0  X  =  0 

(S(x)Y")  =  (S(x)Y")'  +  PY'  =  0     X  =  L 

Stiffness  and  mass  distributions  for  the  column  are 

S(x)  =  E(x)I(x) 
M(x)  =  p(x)A(x) 

where  I(x)  is  the  second  moment  of  area  about  the  neutral  axis. 

It  can  be  shot.;m  that  assuming  similar  cross  sections  provides 
a  relationship  between  I(x)  and  A(x) .   If  f(x)  represents  the  depen- 
dence of  one  dimension  of  the  cross  section  with  the  position  along  the 
beam,  then  for  some  reference  area  A  , 

A(x)  =  Aq  f2(x)  (5.2.2) 

Choosing  a  reference  area  moment  I^  related  to  A   through  a  radius  of 
gyration  k   for  the  cross-sectional  shape  such  that 

it  can   easily   be   shovira    that 

^n 

I(x)    =  ~  a2(x) 
^0 

A  dimensionless  form  of  the  differential  equation  is  obtained 
by  introducing  the  following  dimensionless  quantities 
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t  =  x/L 
n  =  Y/L 


a(x)  =  A(x)/A^ 


e(x)  =  E(x)/Eq 


p(x)  =  R(x)/Rq 


in(x)  =  M(x)/M, 
s(x)  =  S(x)/S 


where  A  ,  E  ,  p   designate  some  reference  value  and 


^0  =  Vo 


and 


m(x)  =  a(x)r(x) 
s(x)  =  a2(x)e(x) 

Design  variables  and  the  given  geometric  constraints  (5.1.1)  are  also 
converted  to  dimensionless  form: 
\   f  a(x)  <  a^ 

e^  <  e(x)  <  e^ 

where  the  transformed  bounds  have  been  divided  by  the  appropriate  refer- 
ence value.   Replacing  variables  in  (5.2.1)  by  their  dimensionless 
equivalent  gives 


106 


_d 
d 


^  (s(t)n)  +  An  +  k  ^  (n  /    m(5)dd=  0 


n  =  n  =  0  t  =  0  (5.2.3) 

(s(t)n)  =  -^  (s(t)n)  +  An  =  0        t  =  1 
at 


where 


A 

PL2 

M^gL^ 

K 

^0^0 

Integrating  by  parts  and  using  the  given  boundary  conditions, 
equivalent  forroulations  to  (5.2.3)  are  obtained  which  are  used  to 
convert  the  governing  equation  to  a  state  component  representation. 
Equivalent  formulations  are: 

1 

(s(t)n)  =  Aln(l)  -  n(t)]  +  k  /       ni(5)[n(S)  -  n(t)]d5 

--   (s(t)n)    =  -   An  -  kn   /        mU)dZ  (5. 2. A) 


ft   ('^"  Tt  ^■"^'^)-'^^)   =  ''^^'^ 


Let 


X  =  n  Deflection 

X2  =  n  Slope 

x„  =  s(t)ii  Moment 

X,  =  n  -T-  (s(t)ri)  Shear  ^  Slope 
4      dt 
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then,  from  this  choice  of  state  components  and  the  last  of  the 
equations  (5.2.4)  it  follows  that 

^1  =  ^2 
^2  "  ^3/s(t) 
x^  =  x^x^ 
X,  =  km(t) 

To  complete  the  transformation  of  the  structural  problem  to  an  optimal 
control  problem,  the  design  parameters  are  treated  as  the  components 
of  control  _u(t)  where 

u^(t)  =  a(x(t)) 

u^Ct)  =  e(x(t)) 

U3(t)  =  r(x(t)) 

Hence  with  the  definition  of  control  _u(t)  and  the  representation  of  the 
system's  state  by  a  four- dimensional  vector,  the  governing  equations  are: 

x-j^  =  x  X,  (0)  =  0 

x^  =  x^/u^u  x^CO)  =  0 

(5.2.5) 

x„  =  x„x,  x„(l)  =  0 

J     2  4  5 

^4  "  ^V3  ""4^^^  "  "^ 

To  complete  the  conversion  of  the  structural  problem  to  an 
optimal  control  problem,  the  specified  x^eight  condition  must  be  con- 
sidered.  In  tei'ms  of  the  specific  mass  distribution  M(x) ,  that  condi- 
tion is 
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L 

W  =  g  /  M(x)dx 

0 


or 


1 

W  =  p  A   gL  /   in(t)dt 
0 

This  introduces  yet  another  dimensionless  parameter  to  the  system  in 

that 

1 
1  =  y  /   u  (t)u  (t)dt 
0 


where 


Pq^qSL 


w 

With  this  last  development,  the  buckling  problem  can  be  stated 
with  more  mathematical  precision.   Find 

UqP^  =  ARGMAX  [-x^(l)]   ,    -x^(l)  =  A 

_u  £  n 

subject  to: 

(i)   differential  constraints  x  =  f(x,u) 


(ii)   kinematic  boundary  conditions     x  (0)  =  0 

x^CO)  =  0 

natural  boundary  conditions       x„(l)  =  0 

x^(l)  =  -> 

(iii)   hard  geometric  constraint         a^  <  u, (t)  <  a 

L  -   1^   -   U 

(iv)   hard  material  constraints         e   <  u  (t)  <  e 

Li       ~  Z  ~      U 

1 
(v)   subsidiary  constraint  p  /   u  u,dt  =  1 

0 
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A  cursory  examination  of  this  problem  statement  reveals  the  inherent 
difficulty  of  optimal  eigenvalue  problems.   If  the  eigenvalue  A  is  the 
cost  functional,  the  cost  functional  appears  in  the  governing  equation 
and/or  the  boundary  conditions.   For  the  state  components  selected, 
the  eigenvalue  X   appears  only  as  a  boundary  condition.   However,  as 
stated  in  (5.2.3),  the  eigenvalue  appears  in  both  the  equation  and 
boundary  conditions. 

5. 3  Analysis  of  the  Problem 

In  the  preceding  section  the  governing  equation  for  column 
buckling  is  expressed  as  either  a  second,  third,  or  fourth  order  equa- 
tion \'n-th   appropriate  boundary  conditions.   State  coipponents  are  then 
defined  from  structural  quantities,  one  of  which  is  the  shear  divided  by 
the  slope.   Since  the  slope  vanishes  at  t  =  0,  it  must  be  sho^'/n  that 
this  state  variable,  x. (t),  is  not  indeterminate  at  the  point  in  ques- 
tion.  Moreover,  in  order  to  use  the  theoretical  techniques  developed 
in  preceding  chapters  the  first  eigenvalue  (cost)  must  be  expressed  as 
some  functional  of  the  state  and  design  parameters. 

According  to  Bolotin  (1963,  p.  22),  conservative  systems  are 

described  by  self-adjoint  boundary  value  problems.   Classical  elastic 

stability  tb.eory  is  conventionall}'  restricted  to  conservative  problems 

in   which  the  buckling  load  is  the  fundamental  eigenvalue.   It  is  shown 

in  what  follows  that  the  column  buckling  problem  is  self-adjoint  and 

that  the  eigenvalue  may  be  obtained  via  the  total  potential  energy  of 

the  system  from  the  Rayleigh  quotient.   Sufficiency  for  x,(t)  to  be  deter- 

4 

iiiinate  at    t  =  0    is    that    the   problem  be   a   self-adjoint   eigenvalue  problem. 
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Consider  a  dependent  variable  in  the  third-order  representa- 
tion of  the  governing  equation  (5.2.4)  and  appropriate  boundary  con- 
dition.  If 

z  =  n 
another  valid  representation  of  the  system  is 


~   (sz)  +  (X+k  /  in(0dC)2  =  0 


''  0 


z  =  0   ,  t  =  0 

s(t)z  =  0   ,  t  =  1 


Using  the  notation  of  Lovitt  (1924),  who  treats  the  general  problem 
(pu')'  +  (q  +  Xr)u  =  0 


0- 


(puu') 
the  corresponding  quantities  are 


u  '^  z 

X  '^  t 

p(x)  '^^  s(t) 

1 
q(:0  '^^  k  /   in(^)dC 

t 

r(x)  ^   1 

UTien  written  as  a  linear  differential  operator  L(  )  and  its  adjoint 
L  (  ),  the  operator  L(  )  is  said  to  be  self-adjoint  if  L  (  )  =  L(  ). 

Tf  the  TPBVP  is  to  be  self-adjoint   for  any  admissible  function 
(f)  and   tjj,  then  by  Green's  Formula 


Ill 


/  {t/'L(<f>)  -  (^L  (4))}  dx  =  B((f>,^) 


=  0 


where  B(4),iIj)  is  the  bilinear  concomitant  of  L(  ).   In  terms  of  a  general 
second  order  equation 


L(*)  =  Pji-t-"  +  P24''  +  (q  +  Ar)<}) 


L  W    =  (p^,|;)"  -  (p^^-)'  +  (q  +  Xr)i, 

The  necessary  and  sufficient  conditions  that  L(  )  be  self-adjoint  is 
that  p„  =  p'.   If  L(  )  is  a  self-adjoint  operator  and  the  bilinear 
concomitant  vanishes,  the  problem  is  said  to  be  self-adjoint.   For  the 
general  operator, 
b 


B(4.,i>) 


=  p,  (cj)'./^  -  (f.^'') 


or  in  terras  of  Lovitt's  equation,  which  obviously  has  a  self-adjoint 
operator. 


B(v^^) 


=  p(<J.'iJ;  -  <}.i|;') 


it  is  seen  that  the  problem  is  self-adjoint  if  the  boundary  conditions 
of  the  adjoint  and  given  systems  are  identical. 

There  is  a  well-developed  theory  associated  with  self-adjoint 
eigenvalue  problems,  where  the  eigenvalue  problem  itself  is  the  result 
of  minimizing  an  energy  functional.  Along  this  line,  Lovitt  shows  that 
under  the  requirements  that  p(x)  >  0  and  is  piecewise  smooth  in  (0,1), 
and  that  q (x)  is  a  well-behaved  function  which  is  finite  everywhere  on 
(0,1),  there  are  an  infinite  number  of  real  positive  eigenvalues  which 
may  be  arranged  as 
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0<A,  <A^<.,.<A   <... 
12  n 

If  p (x)  and  q(x)  are  specified  functions  which  satisfy  the  above  con- 
ditions, the  fundamental  eigenvalue  is  given  by 

A   =  Min  [D(u)] 
u(x) 

1 

D(u)  =  /   {p(x)(u')2  -  q(x)u2}  dx 
0 


where 


u  £  C2 


1 
puu'    =  0 
10 


and 


1 
/   r(x)u2dx  =  1 
0 

2 
Lovitt  shows  for  u  e  C   and  p(x),  q (x) ,  r(x)  bounded,  that  A   is  finite. 

Higher  ordered  eigenvalues  are  obtained  by  using  the  ortho- 
gonality property  of  eigenfunctions   as  additional  subsidiary  condi- 
tions.  For  example, 

A   =   Min   [D(u  )] 
n      ,  ,      n 
u    (x) 
n 

1 

D(u  )  =  /   {p(x)(u')2  -  q(x)u2}  dx 
0         "  ^ 
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for 


u  e  C 
n 


,1 
'0 


pu  u'     =  0 

n  n 


1 

/   r(x)u  dx  =  1 
0 


and 


1 
/   r(x)u  u.dx  =0        1  =  1, . . . ,n-l 
0 

The  immediate  use  of  the  self-adjoint  property  of  the  column 

buckling  problem  is  that  the  fundamental  eigenvalue  is  finite.   This 

is  used  to  prove  that  x,(0)  is  determinate.   Recall  the  definition 

.-1  d  ,  ... 
X,  =  n   — -  (sn) 
4       dt 

From  the  third-order  formulation  of  the  governing  equation  (5.2,4) 

1 


n   J£   (sn)  =  -  A  -  k  /    m(C)d5 


.-1  d 

:  (.sn;  =  -  a  -  tc 

t 

such  that 

1 

x,(0)  =  -  A  -  k  /  m(C)dC  =  -  3 
0 

where  p  is  a  positive  constant  to  be  determined.   By  substituting  the 

specified  weight  subsidiary  condition  for  the  integral 

x^(0)  =  -  (A  +  k/p)  =  -  B 
Both  k  and  p  are  specified,  finite  parameters  of  the  system.   If  A  is 
finite  then  x, (0)  is  also.   In  fact,  the  buckling  load  X    is  the 
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fundamental  eigenvalue  A   which  is  shown  by  Lovitt  to  be  finite  for 

the  conditions  prescribed.   Hence,  ^,(0)  is  determinate. 

In  the  literature  of  classical  elastic  stability  problems, 

the  eigenvalue  to  be  optimized  is  expressed  as  a  Rayleigh  quotient. 

This  quotient  can  be  found  from  the  total  potential  energy.   For  the 

general  equation  of  Lovitt,  the  quotient  is 

1 
/   {p(u')2  -  qu2}  dx 
0 


A  = 


1 

/   ru^dx 
0 

or  in  terms  of  the  column  buckling  parameters  and  components  (5.2.5) 

1  1 

/   {x2/s(t)  -  k  /  m(5)d5  x2}  dt 
0    ^  t  ^ 


A  = 


1 

/   x2  dt 
0 

On  reversing  the  order  of  integralon  for  the  double  integral  contained 
in  the  numerator,  and  introducing  design  parameters,  the  Rayleigh 

quotient  becomes 

1  t 

/   {x2/u2u   -  ku  u   /   x2(5)dU  dt 
_  0__2_  0 


/'  x|  dt 
0 


The  normalization  condition  of  Lovitt  requires  that  the 
quotient's  denominator  equals  unity  for  all  admissible  solutions,  or 

1 

/   x2(t)dt  =  1 
0 
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This  particular  normalization  is  selected  in  order  to  obtain  the 
simplest  expression  for  the  functional  A.   Also,  recall  that 


x^(0)  =  -  (X  +  k/y) 


and 


x^(l)  =  -A 

It  is  seen  from  these  boundary  conditions  and  the  definition 

X.  =  ku  u   that  the  subsidiary  weight  condition  can  be  ^-rritten  as  a 

mixed  boundary  condition  on  x. (t) : 

x^(l)  -  x^(0)  =  k/p 

This,  together  with  converting  the  normalization  condition  to  another 
state  component  with  boundary  conditions  specified  at  both  ends, 
transforms  the  problem  to: 


A..   =  Max  {J[u]} 
Max  '— ' 

u  e  Q 


subject  to 


1 
J[u]  =  /   {x^/u^u-  -  ku  u  X  }  dt 


i,  =  x^  x^(0)  =  0 


i^  =  ^2-''l'^2  ""2^°^  "  ° 


X3  =  x^x^  X3(l)  =  0 


x^  =  x2  x^(0)  =  0 

x^d)  =  1 


116 


where  _u  e  f^  means  that  ueC^    and  satisfies  the  geometric  control 
inequalities.   For  reasons  as  yet  undetermined,  the  quasilineariza- 
tion  algorithm  did  not  converge  when  the  latter  two  state  boundary 
conditions  were  used,  even  though  they  belong  to  the  class  of  problems 
to  Vvfhich  the  algorithm  is  theoretically  applicable.   This  might  be  due 
solely  to  the  mixed  boundary  condition  being  dependent  on  the  value  of 
the  cost  functional  being  minimized. 

On  failing  to  obtain  a  solution  by  minimization  of  the  Rayleigh 
quotient,  the  next  attempt  was  suggested  by  the  state  boundary  condi- 
tions.  From  the  condition 
x^(l)  =  -  A 

and  definition  of  A,  whereby  A  >  0,  it  follows  that 

Min  [x  (1)]  =  Max  [A] 
_u  u. 

Tlius,  the  problem  can  be  cast  as  a  Mayer  type  of  optimal  control  problem: 

Vx=  '""l     f^(^>] 

^^   =  x^  x^(0)  =  0 

X  =  X.-./U-U  x„(0)  -  0 

S  =  V4  X3(l)  =  0 

X,  =  ku  u  X. (0)  =  -B 

4     1  3  4 

where  B  is  a  positive  constant  to  be  determined.   The  quasilineariza- 

tion  algorithm  did  not  converge  for  this  formulation  of  the  problem 

also,  again  for  undetennined  reasons. 
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Consider  from  the  choice  of  state  variables  that 

X,  =  km(t)     and     x, (1)  =  -  A 
4  4 

In  an  earlier  derivation,  3  was  defined  such  that 

1 

X,  (0)  =  -  A  -  k  I   in(t)  dt  =  -  3  <  0 
0 

In  another  form, 

1 
-  6  +  k  /  m(t)  dt  =  -  A 
0 

(Constant)  +  (Weight)  =  -  (Load) 

or  in  terms  of  the  weight 

A  =  3  -  kW 

Therefore,  if  the  weight  W  is  minimized  subject  to  constraints,  this 
is  the  same  as  maximizing  A,  the  load,  subject  to  the  same  constraints. 
Notice  that  constant  3  is  an  undetermined  boundary  condition,  i.e., 
a  free  condition,  so  that  the  corresponding  adjoint  variable's  boundary 
condition  is  known. 

At  this  point  the  results  of  a  doctoral  dissertation — Salinas 
(1968) — were  used  in  another  formulation  of  the  same  problem.   Salinas 
considered  the  application  of  the  energy  method  to  self-adjoint  prob- 
lems.  He  showed  that  the  follov;ing  problems  are  identical: 

(i)   Maximum  load,  specified  weight 
UqP^  =  ARGMAX   {Jp[u]} 

J  =  the  buckling  load  cost  functional 
J  =  a  constant,  the  weight 
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(ii)   Minimum  weight,  specified  load 

UqP^  =  ARGMIN  {J^^  [u]} 
_u  e  fi 

J   =  the  weight  cost  functional 

w 

J   =  a  constant,  the  load 

Both  problems  are  subject  to  the  same  differential  constraints, 
boundary  conditions,  and  control  admissibility  constraints.   A  concise 
way  of  stating  the  equivalence  is 

_ii  ^„  =  ARGI-lAX   {Jp  [u]  :  X^  specified} 
—OPT       ^     P  —     w 

=   ARGMIN   {J   [u]  :  J   specified} 
_u  e  fi 

The  equivalence  of  these  two  problems  is  used  to  state  the 

column  buckling  problem  as 

UqP^  =  ARGMIN  [J^^] 


\  =   ^'/^    ^^^3 

dt 

^1  =  ^2 

x^(0) 

=   0 

^2  =  ^3''"i"2 

X2(0) 

=   0 

^3  =  ^2^ 

X3(l) 

=   0 

h  =  ''"l^S 

x^(l) 

=  - 

A 
s 

where  the  subscript  s  on  the  eigenvalue  indicates  a  specified  load. 

With  reference  to  the  statement  of  the  problem  as  a  Mayer  type,  notice 

that 
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J.J  =1   [x,  (1)  -  X,  (0)]  ->  Min   [x.  (1)] 

u  £  i2 

with  X, (0)  =  -  3,  a  free  boundary  condition.   Furthermore,  since 

X, (1)  =  -  A    and   A  E  J 
observe  that 

!  \  =  ^  -  ^p 

Admissibility  of  control  (ii  e  9.)    requires  that  u.(t)  is  piecewise 
continuous  through  the  second  derivative,  and  that  the  control  inequal- 
ity constraints  are  not  violated.   For  every  value  X   assigned  a  priori 

to  the  eigenvalue  A  there  is  a  solution  (u^    ,  ^         )  that  minimizes 

OPT    OPT 

the  weight  with  respect  to  u  e  f2  and  satisfies  the  eigenvalue  problem. 
However,  an  arbitrary  A   need  not  satisfy  the  specified  weight  condi- 
tion and  in  general  does  not.   To  fulfill  the  given  total  weight  require- 
ment the  value  of  A  was  adjusted  such  that 


J„  [u^j   ]  =  y/  u^u^dt  =  1 


OPT      0 

When  the  cost  is  not  only  a  minimum  but  is  also  equal  to  unity,  all 
conditions  of  the  buckling  problem  as  stated  in  section  5.2  are  satis- 
fied.  By  Salinas'  theorem  the  A   causing  this  satisfaction  is  the 
maximum  lowest  eigenvalue,  i.e.,  the  maximum  buckling  load.   This  effect 
is  illustrated  in  Figure  5.2  and  is  mathematically  stated  by 

Max     r  -^v^^p^      w  "^Qp^ 


where 


LL     =  ARGMIN  {J  [u]  :  A  =  A  } 
OPT    u  e  f2  ^ 
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5 .4  Application  of  the  Maximum  Principle 

In  the  preceding  analysis  difficulties  associated  with  the 
optimal  eigenvalue  problem  are  discussed.   Circumvention  of  such  dif- 
ficulties is  described  for  those  problems  that  are  self-adjoint. 
It  is  shown  that  the  column  buckling  belongs  to  that  class  of  prob- 
lems.  However,  before  the  theoretical  developments  presented  in  the 
third  chapter  may  be  applied,  the  particular  forms  of  parameters  and 
functions  must  be  identified.   Accordingly,  for  the  optimization  of 
the  column  buckling  load 

X  '^.  (x^  x^  X3  x^)^ 

u^  (u^  U2  U3)' 

1 

J  '^'  J,  =  /   L  (x,u)dt 
W    p   U 


Lq  (^I'Ji)   "    ^"1^3 


r 


f(x,u)  = 


"3/"i"2 


^2^4 


^u^U3 


Boundary 
Conditions 


x^CO) 
X3(l) 

1x^(1) 


0 
0 
0 

-A 


122 


(a^-u^) (u^-a^) 


'i'j^Cii.ii)  '^  <  -  ^^^'"^2^ '^"'r^'l?      '  ^£  \  ^2  (5.A.1) 


where  the  maximum  value  of  the  As  is  that  value  such  that 

1 

Jy  "  "  ^   "l^S^t  =  1 
0 

Based  upon  these  particular  forms  for  the  general  parameters, 
the  variational  Hamiltonian  is 


H 


{p^x^  +  V2^^/^i^2   •'  P3^2^  +  Pa'^'^i^s^  ""  ^"1^3 


Adjoint  variables  p^(t)  are  defined  from  the  Hamiltonian  H  which  con- 
tains the  constraint  functions.   Since  the  constraint  functions  (J)   are 
independent  of  x  for  the  case  of  only  geometric  constraints 

p  =  -  h"  =  H 

From  this  result  and  from  the  state  boundary  conditions  it  follows  that 

P-L  =  0  p^(l)  =  0 

P2  =  -Pl^2  -  PS^        P2^^^  ^  ° 

P3  =  -^2^-1-2  P3^°^  =  ° 

P4  "  ""^3^2  Pa^^-*  "  ° 

Variable  p,(t)  may  be  integrated  by  inspection: 

p  (t)  =  0     for  all     t  >  0 

which  immediately  reduces  the  number  of  adjoint  variables  to  be  deter- 
mined and  simplifies  the  Hamiltonian  as  well.   Thus, 
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H  =  -  {p  X2X^  +  p^x  /u^u  +  (kp  -ij)u  u  } 
where 

P2  =  ~  P3^4  V>2^^^    "  ^ 

P4  =  -  P3^2         ^4*^°^  "^  ° 
Optimal  control  as  determined  by  the  PMP  stated  in  (3.1.8)  is 

Uqp^  =  ARGMIN  [H(Xqp^,u,p)]  (5.4.2) 

where  the  admissibility  restriction  in  Pontryagin's  derivation 
requires  only  that  ja(t)  be  piecewise  continuous.   Applying  Salinas' 
theorem  further  restricts  the  control  to  _u  £  C^.   Considering  for  the 
moment  that  no  further  constraints  on  the  control  exist,  the  uncon- 
strained u^-prp  is  that  choice  of  (u,  (t) ,  u„(l),  u„(t))  which  minimizes 
the  variational  Hamiltonian  H.   Recognizing  that  the  physical  quantities 
represented  by  u  have  meaning  only  if  they  are  positive,  assume  this  to 
be  the  case.   With  non-negative  controls,  the  state  equations  and 
boundary  conditions  jointly  imply  that 

x.(t)  >  0     i  =  1,2,3 

x^(t)  <  0 

Furthermore,  by  means  of  (5.4.2)  and  the  expression  for  H 

(p^x^)  >  0  ->  U2  =  0 

(p^x^)  <  0  •>  u^  =  " 
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and 

(kp,-y)  >  0  ->  u  =  °=> 

(kp^-y)  <  0  ->  u^  =  0 

Control  u   remains  to  be  determined. 

Assuming  that  u„,u   are  known  quantities,  then  u^  must  be  chosen 
independently  to  minimize  the  variational  Hamiltonian.   In  that  context 

H^  =  -  {-2p2X^/u3u2  +  (kp^-y)u^}  =  0 

H     =  -  {6p„x_/u';^u„}  >  0 
u  u         2  3   12 

The  unconstrained  optimal  control  variable  u   is  defined  to  satisfy  the 
first  condition;  the  second  condition  indicates  that  p^(t)  must  be  nega- 
tive for  H(x^p^,u^p^,£)  to  be  a  minimum.   So, 

"l  -  V^kp^-p)u2U3J 
and 

PjCt)  <  0 

which  leads  to  the  following  argument  on  the  signs  of  p.(t): 
P2  >  0  ^     p^Ct)  <  0      and   p^CD  =  0 
p^Ct)  >  0  ^     p^  =  -V-^i^     and   x^(t)  <  0 
P^^  <  0  -^     p^  =  -P^x^   and   x^Ct)  >  0 


p^(t)  <  0  ^     Pa  <  0      3nd   Pa^O)  =  0 
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It  is  noted  that  with  these  signs  if  (kp  -y)  <  0  the  optimal  uncon- 
strained control  is 


^PT  =  \  "^     / 


vo  J 


where  the  existence  of  u   is  discussed  in  the  next  section.   What  is 
important  is  that  the  solution  requires  a  material  of  infinite  strength 
and  zero  density.   Some  structural  analysts  commonly  refer  to  this 
hypothetical  material  as  Bolognium, 

In  the  course  of  an  actual  design  process,  control  bounds  corre- 
sponding to  the  constraint  functions  are  specified.   On  the  basis  of 
the  preceding  analysis  with  given  bounds  on  controls,  optimal  u„(t)  and 
u,,  (t)  are  bang- bang  controls  defined  by 

(5.4.3) 


U2(t)  =  %(ej+eu)  +  %(e-^-eu)  SGN  (P2>:3) 
U3(t)  =  %(r^+r^.)  -  hi.^^-^^    SGN  (kp^-y) 
where  for  an  arbitrary  argument  z  of  function  SGN  (  ) 


(5.4.4) 


r+1 

SGN  (z)  =  <   0 
v-1 


z  >  0 
z  =  0 
z  <  0 


Depending  upon  the  magnitude  of  u  '  vis-a-vis  a  and  a  ,  there  are  three 


possibilities  for  u^   at  each  point  (t,2£(t)  ,£^(t))  : 
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(X)   u^  <  a^ 


^PT 


\'' 


"2 


^2'' 


(ii)   a   <  u   <  a      u 

^   ■'    L     1     U     -^PT 


^3  J 


(iii)  a^  <  u 


^^,^ 


i!0?T=\ 


u. 


0^3^ 


That  class  of  problems  having  no  nii:ced  constraints  requires  computing 
the  Lagrangian  multiplier  functions  jj(t)  of  (3.3.4)  only  as  a  check. 
Each  individual  multiplier  function  must  be  non-negative  for  all 
0  <  t  <  1.   A  positive  value  in  a  mathematical  programming  application 
indicates  that  the  usable,  feasible  direction  is  attempting  to  leave 
an  active  constraint  surface.   With  the  similarity  betx/een  this  and 
Pf'lP  demonstrated  in  Chapter  III,  the  existence  of  a  negative  multiplier 
function  indicates  a  similar  inconsistency  in  active  constraints  and 
the  constrained  direction  of  steepest  descent  of  H  with  respect  to  _u. 
In  the  usual  application  of  PMP  to  a  problem,  the  constraints 
are  adjoined  to  the  cost  function  and  another  functional  defined: 


H*  =  -  H  -  p^(t)  l(x,u) 
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Each  constraint  <^gi^,u)   which  is  active  such  that 
'i'^^Cx.u)  =  0 

is  solved  for  a  component  of  control.   The  resulting  non-zero  multi- 
plier function  P„(t)  must  be  determined  from  the  necessary  condition 

* 

H   =0 
u. 

This  \i    (t)  is  used  in  turn  to  determine  the  effect  upon  the  system 
from  the  solution's  lying  on  a  mixed  constraint  surface  v/hich  by 
definition  depends  upon  the  state.   This  effect  is  manifested  by   the 
adjoint  differential  equation  _g^  =  -  H  which  can  be  written  as 

£  =  -  (H^  +  ^) 

If  constraints  <}>   depend  only  upon  the  control,  the  second  term  is 

T 
unnecessary  because  A   -  0. 

Such  is  the  case  of  the  column  buckling  problem  possessing  only 

control  bounds.   For  this  problem,  xjhenever  the  constraints  are  active 

the  multiplier  functions  r^eterm-ined  by  the  above  procedure  are 

-2p  X  /u^u  +  (kp  -vi)u 

1  L     u 

-2p  X  /u2u2 
2  2u2  -  (e^+e^) 

(kp  -m)u 
3     u   L 
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These  same  parameters  will  now  be  determined  from  the  explicit  formula- 
tion of  _p(t)  derived  in  Chapter  III. 

Recall  that  where  I.  denotes  the  set  of  active  constraints  which 

A 

are  r  in  number: 


■"■A  "  ^"l'  °2'*  "  '"'r'^ 


{<{)    ^        ...    <t>    } 
12       r 


K-  f^ii 


^K: 


du. 


1  J 


i  =  1, . . . ,m 

and  that  in  terms  of  ^   and  H,  jj(t)  is  defined  as 

£(t)  =-U,^]       ±^\ 
\    u   ~J  u        — 


(3. 3, A) 


For  the  column  buckling  problem  subject  to  constraints  (5.4.1),  in 
order  to  calculate  jj  the  following  quantities  must  be  defined: 


^1.1 


0 


0 


*l.u=   0    )     ^2,u=    *2,2)     ^3,u=    ' 


0 


0 

^    J 


^Hr,  "^ 


r 


~^^2^3''"l'^2  "*"  ('^P4-1J)"3*^ 

\  =  { Hu2  y  -  { -2P2-3/UH  ^     ' 

0      +  (kp^-y)u^J 


V.^^U3y 


V. 
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with 


*1,1  =  2U-L  -  (a^+a^) 
*3,3  =  2^3  -  (Vr^) 


Consider  first  the  case  where  u   is  not  on  a  control  boundary. 
Only  the  latter  two  constraints  are  active  since  u„  and  u„  are  bang- 
bang  controls,  thus  y,  =  0  and  both  p„  and  p   are  determined  by  (3.3.4) 
From  the  preceding  paragraph's  definitions 


r  =  2   ->  I^  =  {a^,a^}   =  {2,3} 


i>   = 


y  = 


(1  T^  )  = 


u  — 


0  *2,2  0 


0   0 


3,3 


0     0 

0    4) 


3,3_ 


-2,2 


'3,3_ 


u  — 


^2,2 


2,2 


'3,3 


3,3 


"^2,2 

0       ~ 

H 
"1 

0 

^3,3_. 

< 

H 
U2 

H 
L"3j 

H    /(Po   o 

u„  ^2,2 


IV*3,3 
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Hence, 


-P2^^"ii 

^2    2n  -  (e,+e, ) 


^3  =  2U3  -  (r^+r^) 

which  are  identical  to  the  expressions  (5.4.6)  and  (5.4.7)  obtained  by 
the  usual  method. 

Whenever  the  (^^    constraint  is  active,  u^  is  prescribed  and  \i^    is 
nonzero.   If  the  derivation  in  the  preceding  paragraph  is  expanded  to 
include  this  additional  constraint,  the  resulting  expression  for  p,  is 
identical  to  (5.4.5). 

The  result  of  this  analysis  is  a  nonlinear  TPBVP  whose  solution 
is  the  desired  optimal  solution: 


^1   =   ^2 

x^(0) 

=   0 

^2  =   ^3/"i"2 

X2(0) 

=   0 

x^   -  x^x^^ 

X3(l) 

=   0 

^4   =  ^"l"3 

x^(l) 

=   -A 
s 

P2  =  -P3'^4 

P2(l) 

=   0 

P3  =  -P2/"i"2         ^3^°^  "  ° 
P4  =  -P3"2  P4^°)  =  ° 


u^p^  =  ARGMIN  [H(Xqp^,u,2.)] 

u  E  n 


H  =  -  {p  x^x  +  P2^3/^i'^2  "^  (1'^P4~1J)'J2'^3"^ 


\lax  =  ^V^'kiFT^    '■    V^Vt'^    =  ^^ 
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where  u  e  il  requires  that  ii  satisfy  constraints  (5.4.1)  and  Lagrangian 
multiplier  functions  given  by  (5.4.5)  through  (5.4.7)  must  be  non- 
negative.    A  solution  to  this  problem  is  obtained  by  the  method  of 
quasilinearization;  a  subroutine  listing  of  the  equations  required  by 
the  IBM  program  is  given  in  Appendix  C. 

One  additional  preliminary  item  must  be  attended  to  before  pre- 
senting results.   In  order  to  provide  solutions  subject  to  comparison 
standards,  an  eigenf unction  normalization  must  be  prescribed.   Most 
Rayleigh  quotient  problems  in  the  open  literature  normalize  the  eigen- 
function  such  that  the  denominator  of  the  quotient  is  always  equal  to 
unity.   For  the  column  buckling  problem  this  is  equivalent  to 

1 
/  x?(t)  dt  =  1  (5.4.8) 

0 

The  unmodified  results  of  solution  via  quasilinearization  do  not 

satisfy  this  condition;  instead, 

1 
/   x2(t)  dt  =  a2 
0 

where  A  is  some  constant.  A  transformation  which  satisfies  normaliza- 
tion (5.4.8)  such  that  the  differential  equations,  boundary  conditions, 
and  optimal  control  are  unaffected  by  the  transformation,  is  given  by 


Xj^  •*  x^/A 

— 

■^2  ^  ^2^^^ 

P2  -  AP2 

x^  ^  x„/A 

P3  ->  AP3 

^4  "^4 

P4  "     P4 
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With  a  specified  normalization  of  the  eigenfunction,   *"he  state  solu- 
tion 2L-1PT  ™^y  ^^  compared  to  any  other  similar  solution  for  which  the 
normalization  is  kno^^m.   All  subsequent  results  are  normalized  as  above 
to  satisfy  (5.4.8) . 

5.5   Results:   Geometric  Control  Constraints 

Three  dimensionless  parameters  of  the  system  are  derived  in 
section  5.3.   The  first  is   A,  the  eigenvalue  to  be  maximized.   The 
second  is  k,  which  is  a  relative  measure  of  the  weight  of  a  uniform 
column  to  its  stiffness — how  strongly  inclusion  of  the  weight  affects 
the  solution.   The  third  parameter,  p,  is  a  measure  of  the  weight  of 
a  uniform  column  described  by  the  reference  conditions  relative  to  the 
weight  of  the  column  having  distributed  geometry  and  properties. 
Values  of  k  and  y  describe  a  particular  problem  statement  and  must  be 
specified  a  priori. 

Since  the  present  discussion  is  directed  towards  technique,  as 
opposed  to  results  for  particular  geometries,  loading  conditions,  and 
materials,  a  single  value  of  both  p  and  k  is  used  for  all  cases.   Most 
of  the  literature  compares  the  optimal  eigenvalue  to  that  obtained  for 
a  uniform  cylinder  of  identical  v/eight.   This  corresponds  to  a  value 
of  unity  for  p,  selected  so  that  the  results  may  be  compared.   Param- 
eter k  involves  both  material  and  geometric  properties.   For  convenience 
the  material  selected  is  structural  steel  and  the  Euler  buckling  slim- 
ness  ratio  (L/k  )  value  of  60  is  assumed.   With  these  assumptions  the 
values  used  throughout  are 

;i  =  1  k  =  .012 
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At  this  point  a  comment  upon  the  effect  that  the  magnitude  of 
k  has  upon  the  nature  of  uUp„(t)  is  appropriate.  In  the  last  section 
it  was  sho^^m  that  u  and  u  are  bang-bang  controls  dependent  upon  the 
sign  of  switching  functions  (p  x  )  and  (kp  -y) ,  respectively.  Intui- 
tively one  expects  the  strongest,  least  dense  material  to  be  optimal, 
corresponding  to 

(P2X3)  <  0   -V  u^  =  e^ 
(kp^-y)  <  0   ->  u^  =  r^^ 

Also  shown  in  the  last  section  is  that  for  Hamiltonian  H  to  be  minimized 
by  an  unconstrained  value  of  u  =  u  at  some  time,  it  is  necessary  that 
p„(t)  <  0.   Note  that  if  this  does  not  occur  somewhere  in  0  <  t  <  1 
the  control  is  uniform  throughout,  and  the  system  described  is  a  uniform 
column;  unless  A   =  tt^/A  the  solution  of  the  state  is  the  trivial  solu- 
tion.   A  direct  consequence  of  p  being  non-positive  is  that  P/(t) 
is  also.   Hence  any  positive  value  of  k  is  sufficient  to  make  (kp  -y)  <  0, 
so  that  the  ininimum  density  is  required. 

This  insignificant  requirement  that  a  nontrivial,  optimal  solu- 
tion to  the  buckling  problem  exist  provides  an  unusual  result.   l>Jhen 
the  variable  signs  required  for  a  nontrivial  solution  are  used  in  the 
PMP  expressions  (5.4.3)  and  (5.4.4),  the  intuitive  solution  is  shown 
mathematically  to  be  the  optimal  solution.   It  is  for  this  reason  that 
all  subsequent  cases  were  run  with  control  bound  values 


Coding  errors  produced  precisely  this  case,  and  indeed, to 
within  the  numerical  limits  of  the  computer  the  algorithm  converged 
to  the  trivial  solution. 


134 


-u  =  ^-^      "u  =  ^-^ 


e^  =  0.8        r  =  1.0 


Unity  values  were  selected  to  simplify  results;  nonunity  control  bound 
values  were  included  to  reveal  unforeseen  switching  in  U2  and  u^.   None 
were  encountered.   Effectively,  the  problem  as  stated  has  a  single 
control  component  u  (t)  =  a(x/L). 

To  illustrate  the  operation  of  the  quasilinearization  algorithm, 
detailed  results  are  presented  for  one  case.   Data  for  this  case  are: 

(5.5.1) 


a  ->■  CO 
U 


a  =  0.9 

Just  as  in  the  cantilever  beam  example,  the  initial  guess  is  taken 

from  the  uniform  solution  and  must  satisfy  only  the  boundary  conditions. 

Initial  guesses  are  in  terms  of 


such  that 


C(t)  =  cos  (%Trt) 
S(t)  =  sin  (%iTt) 


x^(t)  =  ^  (l-C(t)) 


x^Ct)  =  2S(t) 


x^(t)  =   C(t) 


X, (t)  =  -  X   -  (l-t)k/y 
4         s 


x,(t)  =  -  ^  C(t) 


x^(t)  =is(t) 
x^(t)  =  -  ^  S 


0  <  t  <  1 
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Given  this  initial  guess,  the  value  of  A  which  satisfies  the 

weight  condition  is  A,,   =  2.75.   This  is  demonstrated  by  the  follow- 

Max 

ing  values  of  A   and  J  which  represent  specific  points  on  a  curve 
like  that  of  Figure  5.2: 

s  W 

2.74108         0.998819 

2.74396  0.999210 

2.75000         1.000000 

For  A.,   =  A   such  that  J  ,  =  1  the  convergence  is  measured  in  terms  of 
Max    s  W 

ERROR  as  defined  in  Section  4.4,  and  sho\>m  in  the  following  tabulation: 


ation 

ERROR 

Component 
(Position) 

Cos 

*- 

0 

.696>vi0~^ 

X2(0) 

1.00000 

2.75 

1 

.151-'^10~^ 

X2(.630) 

1.00017 

2.75 

2 

.128^-10"'* 

X2(.555) 

1.00003 

2.75 

3 

.250*10"^^ 

X2(.550) 

1.00000 

2.75 

Deflection  of  the  centerline  is  shown  in  Figure  5.3  which  contains 
three  curves:  the  initial  guess,  a  plot  of  subsequent  iterates  (coin- 
cident for  the  scale  used),  and  the  normalized  solution.   The  optimal 
control  is  displayed  in  Figure  5.4.   Control  functions  for  iterations 
subsequent  to  the  initial  guess  plot  coincident  to  a  single  curve. 
All  non-zero  components  of  the  state  and  adjoint  variables  are  shown 
in  Figure  5.5.   Shear  (represented  by  the  product  x  x.)  is  also  included 
to  give  a  complete  set  of  structural  variables. 
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Iterates 


rmalized   Solution 


0.9 


Figure  5.3   Typical  Deflection  of  Centerline 
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x.(t) 

1 


2. 

^~^^- 
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s: 
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n 
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Figure  5.5   Typical  State  and  Adjoint  Variable  Soltuions 
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Notice  that  for  the  case  where  density  and  modulus  are  equal 
to  unity,  the  condition  on  specified  weight  simplifies  to 

1 

J,,  =  /   u,(t)dt  =  1  (5.5.2) 

W   ^Q   1 

where  implicitly  p  =  1.   As  the  value  of  a   increases,  it  reaches  an 
upper  limit  of  unity:   values  of  a   larger  than  unity  violate  the 

J-i 

ii/eight  condition  (5.5.2).   Thus,  for  no  upper  bound  on  u   there  are  two 
limiting  cases  associated  i-rLth  lower  bound  values, 

a  =  0     and     a  =  1.0 

Values  of  a   not  contained  in  this  range  automatically  violate  some 

given  condition.   Each  value  of  0  <  a   <  1.0  has  a  corresponding  optimal 

eigenvalue  A„   and  design  profile  u,  (t)  '^'  a(x/L)  .   Several  of  these 
°         Max         or  2 

profiles  are  shown  in  Figure  5.6.   In  all  but  the  uniform  column  limit- 
ing case  for  a   =1.0,  the  optimal  profile  has  a  transition  region 


wh 


ere  ^p^  =  u  .   The  point  where  decreasing  u  (t)  intersects  the  lower 


bound  a   occurs  at  a  time  designated  t  .   Optimal  control  is  on  the 

Li  i-i 

lower  bound  for  t  >  t  . 

A  comparative  measure  for  the  effectiveness  of  the  optimiza- 
tion is  the  ratio  of  A.,   to  the  eigenvalue  for  a  uniform,  column  of 

Max 

equal  weight.   This  latter  eigenvalue  is  the  classical  Euler  buckling 

load  denoted  A^,  and  corresponds  to  the  case  a^  =  1.0.   Thus, 
£  L 

A    ^  ^ 

The  dependence  of  the  effectiveness  ratio  A,,  /A„  and  the  lower  bound 
"^  Max  E 

intercept  t   upon  a^  is  shown  in  Figure  5.7.   As  is  typical  of  simple 
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Figure  5.7   Effectiveness  Ratio  and  Lower  Bound  Intercept: 
No  Upper  Bound 
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problems  of  this  type,  the  most  improvement  occurs  for  the  smallest 
excursions  from  the  uniform  case  (.5  <  a   <  1.0).   To  allow  a   to  be 
reduced  from  .6  to  0  provides  much  less  improvement  in  the  cost. 

Observe  that  the  maximum  value  for  a(x)  occurs  at  t  =  0  for 
a   =  0,  and  is  1.333.   Thus  for  a   >  1.333  no  optimal  profile  in 
Figure  5.6  is  constrained  by  the  upper  bound  on  u, (t) .   For  a   =  1.333, 
at  one  point  t  =  0, 

a(0)  =  a^ 
and 

a(t)  <  ay  t  >  0 

Therefore,  instead  of  being  unbounded   in  (5.5.1),  it  is  only  necessary 
that  a   >  1.333.   An  a   tViis  large  is  essentially  unbounded  for  the 
given  weight  constraint. 

Results  similar  to  those  above  have  been  obtained  for  similar 
problems  by  other  methods.   The  next  examples  are  for  a  case  not  yet 
reported  in  the  open  literature.   These  cases  could  not  be  obtained 
without  specifying  both  upper  and  loxjer  control  bounds.   The  first 
example  has  a  lower  bound  specified  with  a  range  of  various  upper 
bounds,  decreasing  to  the  limiting  case  of  a  uniform  solution.   The 
first  example  of  this  case  is 

1.0  <  a   <  1.333 

Since  the  nature  of  the  solution  methods  differ  only  in  details,  no 
solution  methods  are  presented.  Optimal  profiles  are  presented  for 
various  values  of  a   given  a  =  0.   In  Figure  5.8  the  limiting  cases 
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are  shown  with  three  intermediate  cases  as  well.   One  case  is  displayed 
as  a  dashed  line,  since  it  is  obtained  as  an  approximation.   As  a^  ->-1.0 
the  numerical  procedure  becomes  unstable  for  reasons  discussed  in 
Chapter  VII;  the  case  for  a  =  1.05  must  approach  the  uniform  cylinder 
associated  with  a   =  1.0.   Because  the  transition  region  where 
a(x/L)  <  a   is  approaching  a  vertical  straight  line  in  the  limit  as 
a  ->-  1.0,  the  approximation  is  taken  to  be: 

a(x/L)  =  a^  0  <  x/L  <  t^ 

a(x/L)  =a^|l-^^_-^j       t^<x/L<l 

That  value  of  t  which  satisfies  the  specified  weight  condition  under 
the  assumed  form  of  the  approximation  is  given  by 

No  further  use  is  made  of  this  approximation;   its  sole  reason  for 
existence  is  to  provide  an  additional  curve  in  Figure  5.8  to  more  effec- 
tively demonstrate  the  transition  from  the  case  with  no  upper  bound  to 
the  uniform  beam  case. 

The  characteristic  form  of  solutions  is  indicated  by  a  plot  of 
t   versus  a   sho^-m  in  Figure  5.9  together  with  the  effectiveness  ratio 

A,,   /A^.   Just  as  observed  in  the  case  possessing  only  a  lower  bound, 
Max  E 

when  only  an  upper  bound  is  present  the  first  excursions  from  a  uni- 
form column  (1.0  <  a   <  1.2)  provides  the  greatest  cost  improvement. 

In  the  preceding  case,  although  only  t   is  presented,  the 
implication  that  t.  does  not  exist  is  incorrect.   By  observation  of 

Li 
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Figure  5.9   Effectiveness  Ratio  and  Upper  Bound  Intercept; 
Lower  Bound   a^  =0.0 
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Figure  5.8  and  the  definition  of  t   it  is  obvious  that  for  all  cases 

t   =  1.0.   This  occurs  as  a  consequence  of  selecting  a  =  0.0;  in  order 
L  t. 

to  demonstrate  a  more  general  behavior  of  the  solution  process,  a  non- 
zero lower  bound  is  selected  in  order  to  provide  a  non-constant  t  . 

The  third  configuration  of  the  column  buckling  problem  to  be 
considered  is  for  the  control  bounds 

a,  =0.3 

1.0  <  a   <  1.327 

where  the  upper  limit  to  a   is  such  that  for  a   >  1.327  no  upper  bounds 
are  ever  active  as  u^  is  essentially  unbounded  from  above.   Optimal 
profiles  for  the  limiting  and  three  intermediate  cases  are  shown  in 
Figure  5.10.   The  dashed  curve  is  an  approximation  composed  of  three 
straight  line  segments  (numerical  instability  was  again  encountered  as 
the  uniform  column  limiting  case  was  approached) .   Mathematically,  the 
approximation  is 

a(t)  =  a^  0  <  t  <  ty 

(t  -  ty) 

such  that  the  case  which  satisfies  the  specified  weight  condition 
requires  the  control  bound  intercepts  of  the  approximation  to  satisfy 

t.,  +  t^  =  2 
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From  an  extrapolation  of  intercept  values  associated  with  higher  values 

of  a   in  Figure  5.10,  t   is  found  to  be  about  0.985.   The  corresponding 

value  of  t   is  thus  0.882.   Just  as  in  the  preceding  case  it  must  be 

emphasized  that  this  approximation  is  included  only  to  more  clearly 

illustrate  the  transition  of  a(x/L)  from  one  limiting  case  to  the  other 

as  a   decreases  from  1.327  to  1.0. 

Effectiveness  of  design  optimization  is  revealed  in  the  plot 

of  A,,   /A„  versus  a.,  (see  Figure  5.11).   Characteristic  solution  forms 
Max  E         J 

are  also  depicted  there  in  the  plot  of  t   and  t   dependence  upon  a  . 

J_i       u  u 

It  is  to  be  noted  that  t.  deviates  little  from  the  value  of  unity  but 
very  definitely  has  a  functional  dependence  upon  a... 

5 .6   Inequality  Stress  Constraints 

Having  obtained  solutions  for  various  control  constraints,  the 
next  problem  configuration  attempted  involves  a  mixed  constraint. 
Initially  it  was  decided  to  restrict  the  maximum  allowable  normal 
stress  due  to  bending: 

l^B^^^I  ^  ^Max  0  <  X  <  L 

Recognizing  that  for  any  cross  section  the  maximum  value  of  cr  (x)  occurs 
at  the  "outermost  fibres,"  the  distance  d(x)  of  these  fibres  from  the 
bending  axis  must  be  determined.   On  the  basis  of  assuming  similar  cross 
sections,  the  function  f(x)  appearing  in  the  A(x)  expression  (5.2.1) 
also  determines  d(x),  i.e., 

d(x)  =  dQf(x)  =  dQa^(x) 
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Figure  5.11   Effectiveness  Ratio  and  Control  Bound  Intercepts: 
Lower  Bound   a  =0.3 
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where  d   is  some  reference  value.   In  terms  of  Bernoulli-Euler  bending 
equation,  the  inequality  constraint  on  bending  normal  stress  is 


|a,(x)|  = 


Mg(x)d(x) 


I/x) 


< 

-  a  Max 


(EIY")d(x) 


IgCx) 


-  a  Max 


=  |e(x)  d(x)Y"| 


< 

-  a  Max 


By  introducing  dimensionless  variables  and  squaring  to  remove  the 
absolute  value  signs, 


^T  e2(0a(t)n  ;  o^„^^ 


which  is  manipulated  into  the  desired  mixed  constraint  form  when  written 
in  terms  of  state  and  control  variables: 


^H 


°0<  ° 


(5.6.1) 


where 


L  a 


Max 


^O'^O 


Besides  satisfying  the  geometric  constraints,  the  u,  control 
component  must  also  satisfy  a  state  dependent  constraint 

u^(t)  >  (x^/Oq)^^^  >  0 


Results  for  this  case  of  mixed  constraints  are  not  presented.   To  be 
a  valid  demonstration,  the  optimal  control  ought  to  be  constrained  by 
the  mixed  constraint  (5.6.1)  on  at  least  part  of  the  interval  0  <  t  <  1. 
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In  the  present  example,  for  small  a   values  no  solution  exists  that  can 
fulfill  both  the  stress  and  geometric  constraints.   Large  values  of  o 
result  in  all  geometrically  admissible  control  simultaneously  satisfy- 
ing the  mixed  constraint.   Intermediate  values  produce  solution  iter- 
ates in  the  quasilinearization  process  that  were  constrained  by  (5.6.1) 
over  some  interval.   However,  the  final  iterate  u^   (t)  is  not  con- 
strained by  the  mixed  constraint  at  any  point  0  <  t  <  1.   For  this  rea- 
son no  mixed  constraint  examples  are  presented.   Theoretical  consider- 
ations of  other  possible  stress  constraints  are  given  in  the  following 
paragraphs. 

A  more  realistic  stress  constraint  is  to  limit  the  stress 
component  a(x)  normal  to  a  cross  section,  composed  of  a  term  o"  (x)  due 

D 

to  bending  and  a  term  o  (x)  attributed  to  the  end  load  P.   With  refer- 
ence to  the  definition  of  9  shown  in  Figure  5.1,  a(x)  and  its  constitu- 
ents are: 

a(x)  =  OgCx)  +  OpCx) 


Mg(x)d(x) 


^B^^^  =   I^(x) 


,     .     P  cos  0(x) 
Ot,(x)  = 


P'  '      A(x) 

Differential  calculus  and  trigonometry  provide  two  identities  that  are 
used  to  evaluate  cos  (G)  as  a  function  of  x  in  terms  of  the  slope  of 
the  centerline.   They  are 

--  =  tan  9 
dx 

tan^e  +  1  =  sec^g 


^2  1-% 
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From  this, 

Thus,  the  normal  stress  inequality  becomes 

"^"^  =   I,(x)    ^  loo  ['  ^  fe)  j   -  °Max 

By  introducing  dimensionless  state  and  control  variables,  after 
much  manipulation  the  following  mixed  constraint  inequality  is  obtained: 

<||(x,u)  <  0 


Three  parameters  exist  for  this  derivation. 
Max 

-— -  "^   the  Euler  buckling  sliraness  ratio 

—  '^  a  geometrical  slimness  ratxo 
A  simpler  constraint  function  can  be  obtained  by  taking  o  (x)  as  merely 

°p^"^  =ioo 

for  which 


,fe„,  =_(„%,  3,  i(!oj=Ci)„%,,^,„ 
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This  function   can  be  analytically  solved  for  u  whenever 


<|)(x,_u)  =  0.   What  must  be  solved  is  a  cubic  equation  in  u 


1^  _      1, 


(up  3  _  A(u£)  -  B  =  0  (5.6.2) 


-t(^)Vt) 


B  =  X3/aQ 


A  reduced  cubic  equation  with  real  coefficients 

y^  +  py  +  q  =  0 
has  a  discriminant 

A  =  -  108R,     R  =  (^  p)3  +  (~   q)2 


that  determines  the  nature  of  the  roots  to  the  equation.   For, 

A  <  0  ->     one  real,  two  complex  roots 

A  =  0  ->   two  real,  equal  roots 

A  >  0  ->   three  distinct  real  roots 

Consider  a  typical  problem  where  the  cross  section  is  circular,  and 

having  typical  parameter  values: 

d. 


^^0   =    2%              L    = 

30 

°0  = 

.06 

it   follows    from 

1     A    f%\ 

q   =   -   X3/a^ 
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and 


A  =  -  4p3  -  27q3 


% 


(A3/11.664  a^)  x  10   -  x^ 


27   J  ,  c,  „  -,^-5  ,-,   _.? 
a; 


a3  X 10  "/l   -   XI 


Thus  for  typical  parameter  values,  there  is  a  single  real  u,  value 
whenever 


x?   >  a3  X  10   /7  =    4  X  10 


or 

X  >  (a3  xl0"'77)'^  =  6.2  X  10~^ 


In  other  words,  near  the  tip  where  both  t  and  x„  approach  zero,  three 
real  distinct  roots  to  (5.6.2)  exist.   Unless  this  possibility  is 
excluded  by  a  sufficiently  large  minimum  bound  to  u^ ,  the  question  of 
which  root  is  appropriate  must  be  answered.   Having  demonstrated  how 
complicated  the  analysis  becomes  for  this  particular  stress  constraint, 
no  further  considerations  are  presented.   The  developments  are  included 
to  show  that  even  simple  stress  constraints  can  be  too  complicated  for 
ordinary  methods  of  analysis. 


CHAPTER  VI 


FINITE  ELEMENT  METHODS  IN  STRUCTURAL 
OPTIMIZATION:   AN  EXAMPLE 


6 .0  Introduction 

In  order  to  provide  a  contrasting  comparison  to  the  foregoing 
PMP  examples  the  cantilever  beam  problem  is  solved  by  using  finite 
elements  in  conjunction  with  mathematical  programming  techniques.   It 
is  assumed  the  beam  has  a  constant  specified  width  in  order  to  better 
exemplify  the  concepts  contained  in  this  chapter.   Numerical  solutions 
obtained  by  a  feasible  directions  algorithm  are  compared  to  the  PMP 
results.   This  example  serves  to  illustrate  the  principles  used  in  the 
theoretical  derivations. 

6 .1  Finite  Elem.ent  Problem  Statement 

The  same  cantilever  beam  problem  solved  in  Chapter  IV  is  treated 
here  with  different  techniques.   In  what  follows  the  structural  system 
is  modeled  by  and  solved  with  finite  elements.   A  direct  consequence 
of  this  mathematical  representation  is  that  the  optimization  problem  is 
transferred  from  a  vector  space  of  continuous   functions  to  a  space  of 
finite  dimension.   Since  no  contributions  to  finite  element  theory  are 
made,  the  method  itself  is  not  discussed.   The  reader  is  referred  to 
Zienkiewic^  (1971)  for  details  of  the  method.   Bernoulli-Euler  Bending 
Theory  is  applied  to  a  fixed  length  cantilever  beam  of  constant  width 
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and  material  properties,  which  is  symmetric  about  its  central  axis. 
The  problem  is  to  find  the  variable  height  h(x)  that  gives  the  minimum 
tip  deflection  due  to  its  own  height,  subject  to  hard  constraints 
of  maximum  and  minimum  allowable  heights. 

A  higher  ordered  element  is  used  to  obtain  better  accuracy; 
both  the  displacement  and  its  derivative  at  the  nodes  are  prescribed 
to  be  the  independent  variables.   For  this  simple  problem,  n  nodes  divide 
the  beam  into  (n-1)  elements  of  equal  length. 

Hence,  the  displacement  and  its  derivative  are  the  deflection 
and  rotation  of  the  beam's  centerline  at  the  locations  of  the  nodes. 
It  is  also  assumed  that  each  element  is  linearly  tapered,  v^hereby  the 
Bernoulli-Euler  bending  equation  for  element  "i"  takes  the  form 

EI(x)  ^-^  =  M(x) 
dx2 

I(x)  ■-=  j^   wh3(x) 

Position  along  the  beam  ;■:,  height  h(x),  and  displacement  of  the  beam 
y(x) ,  and  slope  Q(x)  of  the  centerline,  are  approximated  by  vectors 
composed  of  discrete  elements: 

X  =  [x^  x^  ...  xj 

h(x)  =  [h(xj  h(x„)  ...  h(x  )f 
• ±     /         n 

T 

i(x)  =  [yCxj)  y(x2)  •••  y(x  )] 
0(x)  =  [e(x.)  o(x„)  ...  e(x  )f 

1     z         n 

T^Tienever  a  height  distribution  is  specified,  a  finite  element  solution 
for  the  beam's  slope  and  deflection  can  be  obtained. 
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6. 2  Mathematical  Programming: 
Gradient  Projection  Method 

With  the  structural  analysis  given  in  the  form  of  a  finite 

element  solution,  the  m.athematical  statement  of  the  optimization  problem 

is 

Given:   (i)   material  properties  p,E 

(ii)   dimensions  L,  w 

(iii)   number  of  nodes  n 

Find:         Mn  [y(x  )  ] 
h(x)     "" 

subject  to  b  <  h(x,)  <  d        i  =  l,...,n 

This  may  be  stated  as  a  nonlinear  mathematical  programming  problem, 

whose  optimality  conditions  are  prescribed  by  the  Kuhn-Tucker  Theorem 

given  in  section  3.2.   For  this  particular  beam  problem  corresponding 

quantities  of  the  theorem  are 

X  '^  h 

F(x)  ^  v(x) 

x=L 

rh(x_j^)  -  d     i  =  l,...,n     j  =■■  i 
g  (x)  -^  { 

*^b-h(x.)      i=l,...,n     j=  i+n 

The  numerical  algoritlim  used  to  obtain  the  optimal  solution 
satisfying  the  Kuhn-Tucker  conditions  is  a  modified  gradient  projec- 
tion method  described  by  Kowalik  ("1970).   Briefly  speaking,  it  uses 
projections  of  the  cost  function  gradient  into  a  subspace  satisfying 
currently  active  constraints.   Starting  from  a  feasible  design  point 
]i  ,  a  sequence  of  feasible  designs  is  generated  such  that 
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h_ ,h,h,...h  ,h   ,,...,h 
-0'-l'-2'    '-q'-q+l'    '-OPTIMAL 

where 


h  ,  ,  =  h  +  a   S 
-q+1    -q     q 

S  =  -  PVF 

a   >  0 

q 

P  =  I  -  N(N^N)~"'"  n'^ 


A  derivation  of  projection  matrix  P  is  found  in  Kowalik  (1970). 
Columns  of  N  are  the  gradient  vectors  of  the  active  constraints;  thus, 
N  is  an  (n  x  r)  matrix  where  r  is  the  number  of  active  constraints  and 
0<r<n.   Ifr=0,  then  by  definition  N  is  a  null  matrix  and  P  reduces 
reduces  to  the  identity  matrix  I.   With  P  so  defined,  those  components 
of   -VF  that  lead  to  a  constraint  violation  are  subtracted  from  the 
direction  of  steepest  descent.   In  this  manner  S  is  the  direction  which 
best  improves  the  cost  function  without  leaving  the  region  of  feasible 
designs  where  g . (x)  <  0.   From  the  form  of  g . ( x)  in  this  particular 
problem,  the  typical  element  N.   of  N  is  defined  as: 

IK. 

"k 
\k  =  ~^T~  i  =  l,...,n     k  =  l,...,r 

i 

where  a  are  those  indices  1  <  j  <  2n  associated  with  the  r  active 
constraints  for  which  g . (x)  =  0.   Hence, 


and     a,  <  n 
k  - 


and     a,  >  n 
k 


N.,     = 
ik 

D 

'^\ 

fl 

1   =a^ 

\.'{ 

ci 

^   =   \ 
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There  is  no  difficulty  in  calculating  the  gradients  required 
by  N;  since  the  g . (x)  are  specified  functions  the  gradients  can  be 
determined  explicitly.   Derivatives  required  to  evaluate  VF  are  not 
so  easily  obtained.   To  circumvent  the  lack  of  an  explicit  formulation, 
the  derivatives  are  obtained  as  follows.   The  i   component  of  VF  is 
found  numerically  by 

gp    F(x  +  Ax)  -  F(x) 


'"i  llAxll 

where    the   only  nonzero   component   of   Ax  is 

1 

Ax.    =  Y7T     X. 
X        10        1 

To  avoid  instability  as  the  solution  is  approached,  the  increment  is 
modified  to 

Ax.  =  ^TT-  kx. 
X   10   X 

where  k  is  the  maximum  percent  change  of  all  components  of  x  from  the 
previous  iteration.  No  numerical  instabilities  were  encountered  v/ith 
this  scheme;  however,  if  too  small  a  value  of  numerical  tolerance  is 
chosen,  the  convergence  criteria  are  never  satisfied  and  the  algorithm 
begins  to  "hunt"  in  the  neighborhood  of  the  optimal  solution.  A  sim- 
plified flowchart  is  shown  in  Figure  6.1. 

To  furnish  a  comparison  standard,  the  Bernoulli-Euler  equation 
for  a  beam  with  a  linearly  tapered  center  section  was  integrated  anal- 
ytically.  This  revealed  some  interesting  qualitative  results  which  are 
discussed  below,  with  quantitative  results  given  in  the  following  sec- 
tion.  First,  the  expression  for  the  second  derivative  of  y(x)  is 
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Figure  6.1   Simplified  Flowchart  of  Algorithm 
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only  piecewise  smooth  due  to  I (x) .   To  integrate  the  equations  for  the 
linearly  tapered  beam  required  dividing  the  range  of  x  into  three 
reqions.   In  the  region  0  <  x  <  x   the  height  h(x)  =  d  and  y(x)  is 
a  quartic  polynomial  in  x.   Variable  height  complicates  the  equation  in 
the  region  x^  <  x  <  x  ;  however,  the  centerline  deflection  is  still 
a  function  of  a  single  variable,  u(x),  defined  as 

u(x)  =  1  -  (f  -  ^)  A 

A  =  2  —  tan  a,    a  constant 
d 


where 


tan  a  = 


2\6J        /X2   x^ 


such  that 


uCx^)  =  ^ 


One  of  the  complicated  terms  of  the  expression  for  y(x)  in  this  region 
of  the  beam  is 

[2  In  (u(x))  -  u(x)] 
which  is  undefined  at  x  =  x„  when  b/d  =  0.   Within  the  third  region 
x„  <  X  <  L,  the  height  h(x)  =  b  and  the  centerline  deflection  is  a 
quartic  polynomial  in  (x  -  x„)  . 

It  is  shown  in  Hornbuckle  et  al.  (1974)  that  in  all  three 
regions  the  deflection  could  be  xjritten  in  nondimensionalized  form 
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^0  "  1^  ^  (if  (^-^-W 

This  result  supports  Dixon's  statement  that  d/L  is  not  a  parameter  of 
the  dimensionless  system.   These  equations  indicate  that  dimensionless 
results  indeed  depend  only  upon  b/d.   However,  for  a  particular  beam  of 
given  dimensions,  the  actual  tip  deflection  does  increase  as  the  square 
of  the  slimness  ratio  L/d. 

6 .3   Results 

Before  any  attempts  were  made  to  obtain  results,  certain  basic 
questions  had  to  be  answered.   First,  the  integration  inherent  to  the 
finite  element  solution  as  described  in  Zienkiewicz  (1971)   is  done 
numerically,  thereby  requiring  determination  of  the  proper  step  size. 
It  was  found  that  with  the  range  of  parameters  used  in  this  problem, 
ten  integration  intervals  per  element  gave  the  desired  accuracy  without 
excessive  computation.   Secondly,  there  was  a  question  of  how  many  nodes 
are  required  to  give  a  valid  representation  of  the  beam.   Because  the 
system  is  so  simple,  as  few  as  tx^o  finite  elements  produce  good  struc- 
tural results  vis-a-vis   tip  deflection.   However,  unless  eight  or  six- 
teen elements  were  used,  the  approximation  of  h(x)  by  h(^)  was  not 
acceptable.   The  results  for  eight  or  sixteen  elements  were  almost  iden- 
tical.  (For  the  latter  two  cases  the  execution  times  on  an  IBM-370/165 
were  15  seconds  and  110  seconds,  respectively.) 
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Basic  operation  of  the  algorithm  and  convergence  to  the  optimal 
solution  is  illustrated  in  Figure  6.2  for  a  typical  value  of  b/d.    The 
profile  is  displayed  at  each  stage  in  the  design  iteration  process, 
including  the  final  design.   A  "stage"  in  the  design  process  occurs 
whenever  another  constraint  is  encountered,  requiring  the  recomputation 
of  VF,  P,  and  S.   To  simplify  the  figure  only  nine  nodes  were  used;  the 
optimal  profile  for  the  beam  having  linear  taper  is  shown  for  comparison. 
The  listing  below  indicates  how  constraints  are  encountered  as  the  numer- 
ical process  approaches  the  optimal  design.   Tip  deflections  for  each 
design  iterate  are  included.   Thus,  for  the  design  process  depicted  by 
Figure  6.2, 

No.  of  Active 
Constraints      Tip  Deflection 

0  .108879      ...   Initial  Design 

1  .027556 

2  .020673 

3  .017103 

4  .016646 

5  .015351      ...   Final  Design 
Comparison  .015332      ...   Linear  Taper  Design 

In  the  earlier  studies,  the  problem  was  found  to  have  a  single 
parameter  b/d;  however,  two  different  beam  profiles  h(x)  were  obtained 
for  this  problem  by  various  authors.   It  vzas  found  with  the  finite 
element  method  that  for  numerical  precision  corresponding  to  slide 
rule  accuracy  the  optimal  profile  is  an  almost  linearly  tapered  beam. 
Greater  numerical  accuracy  provides  a  distinct  curvature  in  the  tapered 
portion  of  the  beam.   This  is  shown  in  Figure  6.3  by  the  two  curves 
generated  with  the  finite  element  method.   The  profile  depicted  by  the 
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dots  is  taken  from  the  work  by  Boykin  and  Sierakowski  (1972) ,  and  is 
presented  for  corroboration.   For  this  latter  curve  the  absolute  magni- 
tude of  the  error  in  satisfying  the  differential  equations  is  less 
than  10 

Since  tip  deflection  is  the  cost  function,  the  tip  deflection 
from  the  finite  element  algorithm  was  compared  to  that  from  the  linearly 
tapered  beam  (Figure  6.4).   Differences  were  less  than  .5%  except  near 
b/d  =  0.   This  is  attributed  to  the  natural  logarithm  term  in  the  linear 
taper  solution  discussed  previously.   No  computer  results  were  obtained 
for  the  linear  taper  solution  with  b/d  almost  zero  because  the  natural 
logarithm  term  became  excessively  large.   This  may  explain  why  all 
earlier  studies  using  this  example  problem  contain  no  results  for 
b/d  <  0.1.   However,  the  finite  element  algorithm  converged  to  a  solu- 
tion for  b/d  =  0  with  no  difficulty. 

Further  analysis  centered  on  a  statement  in  which  Bellamy 
and  West  (1969),   claim:   ".  .  .as  b/d  increases,  the  midsection  of 
the  beam  profile  reduces  in  length  and  increases  in  slope."   Having 
found  quite  close  agreement  with  the  linear  tapared  tip  deflection, 
we  expected  at  least  qualitative  agreement  for  the  angle  of  taper  a. 
Angle  of  taper  for  the  linearly  tapered  beam  and  the  finite  element 
results  are  shown  in  Figure  6.5.   Notice  that  both  curves  approach  zero 
in  the  limit  as  b/d  approaches  one.   Of  greater  importance  is  the  fact 
that  the  angle  of  slope  for  both  methods  is  reasonably  close  only  for 
b/d  <  0.5.   Figure  6.6  shows  that  for  b/d  >  0.5,  if  the  beam  is  divided 
into  sixteen  finite  elements,  the  tapered  section  contains  only  three 
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nodes  because  it  spans  such  a  short  section  of  beam.   This  is  believed 
to  be  insufficient  for  accurate  representation.   The  same  shortcoming 
is  also  apparent  in  Figure  6.7,  which  shows  the  optimal  profile  dimen- 
sions reported  in  the  various  studies.   For  b/d  >  0.5,  the  value  of  r, 

d 

obtained  by  the  finite  element  method  diverges  from  the  other  solution 
method  results. 

To  verify  independence  of  the  solution  to  L/d,  a  series  of  runs 
were  made  with  fixed  b/d  and  the  slimness  ratio  L/d  ranging  from  five 
to  eighty.   Resulting  profiles  at  first  suggested  that  L/d  is  a  param- 
eter, but  closer  analysis  revealed  there  was  an  oversight  in  the  use  of 
numerical  tolerances.   Initially  a  single  tolerance  was  applied  to  both 
geometrical  constraints  and  the  Kuhn- Tucker  test.   The  former  is  inde- 
pendent of  relative  beam  shape,  but  the  latter  is  not;  this  is  implied 
by  the  form  of  the  nondimensionalization  coefficient  given  by  (6.2.1). 

Basically,  the  projection  matrix  P  deletes  those  components 
from  -VF  which  are  normal  to  an  active  geometrical  constraint  boundary. 
Since  all  numerical  schemes  contain  inherent  accuracy  limits,  the  Kuhn- 
Tucker  test  is  considered  to  be  satisfied  when  at  those  points  not  on 
a  boundary 


^y(x^ 

ax. 

1 


"  'g 


Yet  as  the  slimness  ratio  L/d  increases,  for  any  point  on  the  beam  the 
magnitude  of  the  component  of  the  gradient  also  increases.   By  requir- 
ing the  Kuhn-Tucker  test  to  satisfy  the  same  tolerance  for  all  values 
of  L/d,  as  L/d  increases  the  solution  is  effectively  forced  to  become 
more  accurate.   In  fact,  more  accuracy  than  is  possible  from  the 
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algorithm  was  required  and  the  system  began  hunting  around  the  solution, 
failing  to  converge.   By  numerical  experimentation  and  curve  fitting 
it  was  found — Hornbuckle  et  al.  (1974) — that  uniform  accuracy  could  be 
obtained  for  this  particular  problem  when  the  tolerance  used  in  the 
Kuhn-Tucker  test  is 

^0 

^-    . 

N  =  2.77825 

k  =  1.184045  X  10~ 

Subscript  zero  indicates  values  associated  with  reference  length  corre- 
sponding to  Lf-,/d;  the  value  of  five  was  chosen  because  it  represents  an 
approximate  lower  bound  for  the  validity  of  Bernoulli-Euler  bending 
theory.   In  uniform  beams  whenever  L/d  <  5,  the  effect  of  shear  is  no 
longer  negligible.   As  an  example  of  why  care  should  be  exercised  in 
selecting  the  range  of  parameters  to  generate  data,  Dixon  (1967,  1968) 
obtained  results  for  0.2  <  b/d  <  4;  having  assumed  the  effect  of  shear 
negligible,  he  considered  only  the  range  of  L/d  where  it  is  significant. 


CHAPTER  VII 


COMMENTS  ON  NUMERICAL  INSTABILITY  IN 
THE  OUASILINEARIZATION  ALGORITHM 


7.0  Introduction 

Quasilinearization  is  an  indirect  type  of  numerical  scheme  in 
which  a  solution  is  obtained  for  the  differential  equations  that  must 
be  satisfied  by  the  optimal  solution.   Since  a  standard  IBM  SHARE  pro- 
gram is  used  without  modification,  the  method  is  not  discussed.   Inter- 
ested readers  are  referred  to  any  standard  textbook  on  numerical  tech- 
niques, or  to  Ghosh  (1973). 

This  type  of  algorithm  requires  an  initial  guess  of  the  state 
and  then  iterates  from  that  guess  to  the  solution.   Control  functions 
associated  with  the  final  iterate  are  the  optimal  control  functions 
actually  sought.   If  the  process  converges  at  all,  the  convergence  is 
rapid  and  typically  requires  fewer  than  ten  iterates.   However,  the 
paramount  difficulty  of  using  this  method  is  to  select  an  initial  guess 
which  is  "good"  enough  to  yield  convergence.   This  chapter  contains  some 
numerical  data  related  to  specific  cases  of  the  example  problems  where 
convergence  difficulties  are  encountered. 

7.1  Computer  Program  Convergence  Features 

As  explained  in  section  4.4,  convergence  is  measured  by  whether 
or  not  the  iterate  satisfies  the  differential  equation.   In  terms  of 
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the  quantity  ERROR  defined  by  (4.4.2),  in  subroutine  QUASI  the  test  for 
numerical  instability  involves  three  program  parameters: 

ElMAX     the  value  of  ERROR  for  the  initial  guess 

E2MAX     the  value  of  ERROR  for  some  subsequent  iteration 

BLO       a  specified  instability  parameter 

t'Hienever  the  ratio  of  E2MAX  to  ElMAX  exceeds  BLO  the  program  terminates 
on  an  error  exit.   Typical  values  of  BLO  are  10  to  20.   The  primary  mode 
of  convergence  is  for  E2MAX,  which  decreases  with  each  iteration  under 
normal  circumstances,  to  become  smaller  than  a  specified  value  (assumed 
to  be  0.5  X  10   ).   Other  special  modes  of  convergence  contained  in 
QUASI  were  never  used  by  any  of  the  cases  run  for  either  example  problem. 

The  purpose  of  this  test  is  twofold.   In  some  instances  the 
initial  guess  is  not  "good"  enough  and  the  first  few  iterations  possess 
E2MAX  values  greater  than  ElMAX,  but  is  sufficiently  good  for  the  algo- 
rithm to  ultimately  converge.   The  other  purpose  is  to  prevent  excessive 
computation  for  those  cases  where  the  initial  guess  is  so  bad  that  the 
numerical  process  is  divergent. 

Another  related  convergence  test  is  required  by  the  subroutine 
RECIP  which  performs  required  matrix  inversion.   This  test  involves 
a  program  parameter  SUM  which  is  defined  in  terms  of  a  matrix  C   to  be 
inverted,  and  the  inverse  _C   obtained  numerically  by  RECIP.   If  C 
and  C    denote  appropriate  elements  respectively,  then 

SUM  -  E    E   I CT,  C,  .  I 

k    lytj  -^ 
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If  the  matrix  has  been  Inverted  exactly  the  value  of  SUM  is  zero. 
Whenever  SUM  is  less  than  program  parameter  OFFDG  the  inversion  is  con- 
sidered to  be  sufficiently  accurate.   The  value  selected  for  all  cases 

-15 
for  both  example  problems  is  5  x  10    .   Whenever  the  inverted  matrix 

fails  to  satisfy  the  test  on  SUM,  the  subprogram  iterates  upon  the 

inverse  to  improve  it.   If  the  improved  inverse  satisfies  the  accuracy 

requirement,  the  result  is  returned  to  QUASI.   After  five  unsuccessful 

improvement  attempts,  results  of  the  fifth  attempt  are  returned  to  QUASI 

A  singular  matrix  results  in  immediate  termination  of  all  calculations 

via  an  error  exit. 

During  the  operation  of  the  quasilinearization  algorithm  two 

different  matrices  must  be  inverted.   One  of  these,  D. ,  exists  for  every 

element  i  and  each  separate  D.  must  be  inverted.   From  all  of  the  D. 

—1  ~i 

another  matrix  C^  is  defined  which  also  requires  inversion.   The  specific 
form  of  these  matrices  is  not  important;  however,  their  definitions  may 
be  found  in  Ghosh  (1973,  pp.  52-58).   It  is  noteworthy  that  in  the  cases 
run  for  the  example  problems  no  difficulties  were  experienced  in  the 
inversion  of  any  I),  matrix.   All  matrix  inversion  problems  were  asso- 
ciated with  the  C^   matrix. 

7 . 2  Numerical  Instabilities  for 
Cantilever  Beam  Example 

An  initial  guess  for  the  beam  problem  was  obtained  from  the 

solution  for  a  uniform  beam,  that  is,  for  parameter  values  a/c  =  b/d  =  1, 

Cases  were  then  run  for  a  range  on  both  parameters  from  0.1  to  1.   As 

would  be  expected,  no  convergence  difficulties  were  encountered  with 
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combinations  of  parameter  values  corresponding  to  a  nearly  uniform 
beam.   In  fact,  with  the  initial  guess  based  upon  the  uniform  beam, 
only  three  cases  failed  to  converge.   For  this  initial  guess  given  by 
(4. A. 3)  the  following  three  cases  exhibited  numerical  instability. 

a/c  =  0.1         b/d  =  0.3 

=  0.2  (7.2.1) 

=  0.1 
All  other  combinations  of  parameter  values  resulted  in  convergence  by 
the  quasilinearization  algorithm. 

Since  all  of  those  cases  which  did  converge  showed  very  little 
change  in  structural  state  components  x. (t) ,   i  =  1,...,4,  the  uniform 
beam  initial  guess  was  replaced  by  a  guess  based  upon  the  case  where 
a/c  =  0.1  and  b/d  =  0.4.   It  was  suspected  that  the  adjoint  variable 
guesses  v;erG  not  sufficiently  "good."   Hence,  the  equations  (4.4.3) 
were  replaced  by: 

x^(t)  =    .025  *   t^ 
Y.^(t.)    =    .040  *  t 


0  <  t  <  1  (7.2.2) 


x^dt)    =    .100  ■'<    (1-t) 

x,(t)  =  -.400  *    (1-t)^ 
4 

p^(t)  =   5.00   ^"^  t-^ 
p^(t)  =  -1.25   *  t^ 

VJith  this  initial  guess  only  one  case  corresponding  to  the  three  sets 
of  parameter  values  (7.2.1)  failed  due  to  numerical  instability;  the 
first  two  cases  of  these  cases  converged. 
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Having  succeeded  in  forcing  two  additional  cases  to  converge 
by  improving  the  initial  guess,  this  same  approach  was  tried  with  the 
third  case  as  well.   The  initial  guess  was  based  upon  the  converged 
solution  for  the  case  i-d-th  a/c  =  0.1  and  b/d  =  0.2.   Equations  (7.2.2) 
were  replaced  by 

x^(t)  =    .020  *  t^ 

x^Ct)  =    .040  *  t 

x^Ct)  =    .080  *  (1-t)"^ 

.  (t)  =  -.350^-  (l-t)3      °^^^^  (^•2-3) 

4 

p^(t)  =  13.0   A  t^ 

p,(t)  =  -2.00  *   t^ 
4 

With  this  initial  guess  the  case  for  a/c  =  .1   and  b/d  =  .1   converged 

also.   This  last  case  gave  a  complete  data  set  for  parameter  values 

.1  <  a/c  <  1 

.1  <  b/d  <  1 

in  increments  of  .1.   No  cases  for  parameter  values  equal  to  zero  were 
run  since  the  parameters  appear  in  the  denominator  of  various  e:<pressions. 

As  mentioned  earlier,  indirect  methods  such  as  quasilineariza- 
tion  generally  converge  rapidly  if  the  initial  guess  is  sufficiently 
good.   To  try  to  determine  why  the  three  cases  failed  to  converge,  all 
of  the  cases  were  closely  examined.   Two  characteristics  xjere  observed: 
(i)   There  is  little  change  in  the  maximum  value  of  each  separate 
state  component,  regardless  of  the  parameter  values.   Curva- 
ture of  X  (t)  and  x,(t)  did  change  for  lower  parameter  values; 
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this  is  evident  in  different  expressions  used  for  the 
initial  guess  of  these  two  variables.   However,  the  var- 
iable that  actually  caused  the  instability  termination 
was  Po(t)  and  not  a  state  component, 
(ii)   The  maximum  value  of  both  adjoint  variables  increases  with 
decreasing  parameter  values,  with  an  associated  increase 
in  curvature. 

With  p^Ct)  isolated  as  the  source  of  convergence  difficulties  for  cases 
with  small  parameter  values,  Po(t)  was  plotted  for  several  combinations 
of  small  parameter  values. 

The  increased  curvature  mentioned  above  is  shown  in  Figure  7.1. 
All  cases  that  failed  to  converge  exhibit  a  "numerical  instability 
termination"  for  which  p^Ct)  becomes  divergent.   Maximum  error  in  satis- 
fying the  p   equation  occurs  in  the  region  .755  <  t  <  .965  for  these 
divergent  runs.   From  the  figure  it  can  be  seen  that  this  is  the  region 
of  greatest  change  in  the  curvature  of  the  solution  p_(t) .   This  sug- 
gests that  the  linearization  gradients  in  the  quasilinearization  algo- 
rithm may  not  be  valid.   If  that  is  the  case,  then  the  linear  expansion 
about  the  known  iterate  cannot  be  expected  to  result  in  an  improvement. 
A  modification  developed  by  Ghosh  (1973)  to  inhibit  instabilities  asso- 
ciated with  too- large  improvements  has  no  effect  upon  these  numerically 
divergent  cases.   Removal  of  this  possible  source  of  numerical  instabil- 
ity further  implies  the  difficulty  is  related  to  the  linearity  assumptions. 

An  explanation  for  the  increasing  slope  is  available  from  the 
differential  equation  for  the  adjoint  variables,  where 
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Figure  7.1   Solutions  for  Adjoint  Variable  p„(t) 
for  Certain  Cases  of  Interest 
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p   =  (l-t)/u  3u  p  (0)  =  0 

C7.2.4) 

P4  =  -P3  ^4^°^  =  ° 

Since  both  u^  and  u   are  equal  to  their  minimum  allowable  values  near 
t  =  1,  p   increases  with  decreasing  parameter  values  in  this  region. 
However,  the  term  (1-t)  forces  Pn  to  zero  at  t  =  1.   For  unknown  rea- 
sons, the  large  curvature  associated  with  this  effect  seemed  to  pose 
no  difficulties  for  the  numerical  algorithm.   As  seen  in  equations 
(7.2.4),  if  p~(t)  can  be  determined,  then  P,(t)  is  readily  obtained. 
Curves  for  p,(t),  corresponding  to  cases  presented  in  Figure  7.1,  are 
given  in  Figure  7.2.   All  six  cases  exhibit  positive  curvature  and 
possess  none  of  the  inflection  points  seen  in  their  counterparts  of 
Figure  7.1.   The  increase  in  amplitude  of  these  t\sro  variables  for 
decreasing  values  of  b/d  is  illustrated  in  Figure  7.3.   In  this  figure 
it  is  seen  that  with  respect  to  increasing  amplitude  as  b/d  ^  0,  p^(t) 
is  the  critical  variable. 

The  major  disadvantage  of  the  indirect  solution  method  is  that 
little  is  known  about  the  connection  between  convergence  and  require- 
ments upon  the  initial  guess.   For  this  specific  example,  only  one  of 
the  seven  variables  causes  convergence  difficulties,  and  the  result- 
ing numerical  instability  occurs  in  a  region  of  large,  rapidly  changing 
curvature.   No  convergence  problems  are  encountered  when  the  curvature 
of  the  initial  guess  more  closely  resembles  the  final  solution.   All  of 
this  suggests  that  the  question  of  whether  or  not  an  initial  guess 
allows  the  quasilinearization  algorithm  to  converge,  may  depend  upon  the 
respective  curvatures  of  the  initial  guess  and  the  solution. 
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P^(t) 
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Figure  7.2   Solutions  for  Adjoint  Variable  p,(t) 
for  Certain  Cases  of  Interest 
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7. 3  Numerical  Instabilities  for 
Column  Buckling  Example 

Another  type  of  convergence  failure  was  experienced  x^ith  the 
column  buckling  problem.   In  these  cases  which  did  not  converge,  each 
of  the  variables  was  divergent,  increasing  with  each  iteration  begin- 
ning with  the  initial  guess.   All  of  these  numerically  unstable  cases 
were  for  parameter  values  which  resulted  in  a  solution  that  approached 
the  uniform  column  limiting  solution.   Attempts  to  force  convergence 
were  unsuccessful. 

The  first  attempt  was  to  improve  the  initial  guess  used  as  a 
starting  point  by  the  algorithm.   Although  this  technique  succeeded 
with  the  beam  problem,  it  failed  with  the  column  problem.   The  next 
attempt  to  obtain  convergence  v/as  to  employ  the  instability  inhibiting 
modification  developed  by  Ghosh.   The  only  result  was  to  prolong  the 
numerical  divergence.   At  this  point  it  was  suspected  that  perhaps  the 
source  of  difficulty  v/as  related  to  normalization  of  the  eigenf unction. 
However,  when  each  iterate  was  normalized  as  indicated  in  section  5.4, 
the  sole  result  was  that  the  state  variables  maintained  reasonable 
amplitudes  while  the  amplitudes  of  the  adjoint  variables  increased  more 
than  without  normalization.   Doubling  the  number  of  intervals  in  the 
finite  difference  approximation  (from  100  to  200)  was  also  ineffective. 

A  possible  source  of  the  numerical  difficulty  was  the  matrix 
inversion.   For  those  cases  failing  to  converge  each  of  the  jD.  matrices 
was  inverted  accurately  but  the  inversion  of  C^  was  questionable,  as 
indicated  by  the  value  of  SUM.   Since  each  iterate  depends  directly 
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upon  the  inverse  of  C^,  a  degradation  in  C        would  cause  a  correspond- 
ing degradation  in  the  "solution"  at  each  iteration. 

To  illustrate  what  happens  in  these  cases  consider  a  specific 
example:   a^  =  0,  a  =  1.1.   A  plot  of  values  of  SUM  versus  t  ,  the 
intercept  of  the  upper  bound,  is  shown  in  Figure  7.4.   These  data  points 
represent  the  C^   accuracy  at  different  iterations  for  various  initial 
guesses — acceptable  accuracy  is  indicated  by  the  straight  line.   Basi- 
cally, if  an  initial  guess  yielded  a  t   of  approximately  one-half,  each 
iteration  gives  an  increased  value  of  t  with  a  correspondingly  less 
accurate  _C   .   As  t^  approaches  one,  the  solution  approaches  the  limit- 
ing  uniform  column. 

The  result  of  this  is  that  C   may  become  ill-conditioned  if  the 
parameters  are  chosen  such  that  the  solution  is  similar  to  a  uniform 
column.   There  are  no  e:<plicit  data  or  derivations  to  indicate  this  is 
the  situation,  only  the  implicit  suggestion  of  the  data  in  Figure  7.4. 
However,  as  stated  in  the  preceding  section,  there  is  little  known 
about  why  certain  initial  guesses  fail  to  converge  in  the  indirect 
solution  algorithms. 
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CHAPTER  VIII 


CONCLUSIONS  AND  RECOMMENDATIONS 


8.0   Summary  and  Conclusions 

This  dissertation  treats  the  optimal  design  of  elastic  structures 
subject  to  both  hard  inequality  constraints  and  subsidiary  conditions. 
It  is  shown  that  structural  problems  can  be  classified  into  two  general 
types  depending  upon  whether  or  not  the  cost  functional  is  an  eigenvalue. 
For  the  conservative  systems  considered,  which  are  described  by  self- 
adjoint  differential  equations,  it  is  shoism  that  the  minimum  weight 
problem  is  identical  to  the  maximum  buckling  load  problem. 

Pontryagin's  Maximum  Principle  is  analyzed  as  a  nonlinear 
programming  problem  in  Chapter  III.   Specifically,  the  theory  of  the 
gradient  projection  method  is  applied  to  the  maximum  principle,  with 
the  control  subject  to  inequality  constraints.   From  this  an  explicit 
formulation  is  obtained  for  the  Lagrangian  multipliers  that  adjoin  the 
constraints  to  the  variational  Hamiltonian. 

Additional  characteristics  of  the  maximum  principle  result  from 
the  gradient  projection  method  analogy.   In  the  latter,  the  projected 
gradient  is  that  specific  direction  which  minimizes  the  cost  without 
violating  the  constraints  in  effect  at  any  point  in  question.   As 
applied  to  the  maximum  principle,  the  projected  gradient  is 
H  =  -  (H  +  jj  ^u   which  in  u^-space  denotes  the  direction  of  constrained 
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maximum  descent  whenever  jj  is  given  by  equation  (3.3.4).   While  this 
attribute  is  not  pursued  further  in  the  dissertation,  it  may  prove  use- 
ful in  subsequent  applications  of  the  maximum  principle.   To  further 
illustrate  the  nonlinear  programming  principles  used  in  the  theoret- 
ical development,  the  example  of  Chapter  IV  is  approximated  by  finite 
elements,  and  solved  by  a  gradient  projection  algorithm.   Operation  of 
the  algorithm  is  illustrated  and  the  results  compared  to  those  obtained 
with  the  maximum  principle.   When  a  sufficient  number  of  finite  ele- 
ments are  used  to  adequately  describe  the  structure,  the  results  of  the 
two  methods  are  identical. 

A  problem  concerning  the  minimum  deflection  of  a  beam  is  opti- 
mized by  using  classical  maximum  principle  techniques.   It  is  shown 
that  no  finite  solution  is  obtained  for  unconstrained  control;  this  also 
holds  for  the  column  buckling  problem  treated  in  Chapter  V.   Further- 
more, it  is  shown  that  for  the  conservative  system  corresponding  to  the 
buckling  problem,  the  minimum  weight  and  maximum  buckling  problems  are 
identical.   Another  result  of  using  the  maximum  principle  is  that  for 
the  problem  described,  in  order  to  have  a  nontrivial  solution  the  modu- 
lus and  density  must  be  the  maximum  and  minimum  allowable  constant 
values,  respectively.   This  is  a  mathematical  verification  of  an  intui- 
tive result. 

Numerical  instabilities  experienced  by  the  quasilinearization 
algorithm  are  detailed  in  Chapter  VII.   While  no  concrete  results  are 
obtained,  from  the  data  presented   it  appears  that  the  question  of  con- 
vergence is  related  to  the  respective  curvatures  of  the  initial  guess 
and  the  solution. 
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8. 1   Recommendations 

Further  work  suggested  in  the  course  of  dissertation  research 
falls  into  three  categories:   structural  considerations,  mathematical 
theory,  and  numerical  solution  techniques.   Under  the  first  category 
the  most  obvious  extension  is  to  not  invoke  the  linear  bending  assump- 
tion.  Since  the  quasilinearization  algorithm  handles  nonlinear  TPBVP, 
theoretically  there  is  no  need  for  the  assumption.   In  contrast  to  the 
first  of  equations  (5.2.4)  the  complete,  nonlinear  bending  problem  is 
described  by 

(s(t)n)  =  M(t)[l  +  n^]^''^        »  0  <  t  <  1 

1 
M(t)  =  A[n(l)  -  n(t)]  +  k  /  m(5)[n(0-n(t)]d5 

t 

n(0)  =  n(0)  =0 

Another  interesting  possibility  is  to  include  the  effect  of  shear. 

These  structural  aspects  pertain  to  specific  problems.   A  more 
general  result,  and  hence  more  interesting,  is  related  to  expanding  the 
class  of  problems  that  can  be  treated  by  energy  techniques.   It  is  sug- 
gested that  the  mutual  potential  energy  method  and  the  application  of 
adjoint  systems  to  nonconservative  problems  be  investigated  to  see  if 
a  "general  energy  method"  exists.   Such  a  development  would  expand  the 
class  of  problems  now  considered  to  be  amenable  to  the  energy  method 
and  its  well-developed  theory. 

In  terms  of  mathematical  theory,  the  double  optimization 
problem  associated  with  extremization  of  eigenvalues  for  self-adjoint 
systems  is  an  inviting  prospect.   No  significant  contributions  have  been 
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published  in  recent  years.   Any  new  discovery  of  a  general  nature  would 
be  a  significant  development.   Another  useful  area  of  mathematical 
research  would  be  the  further  study  of  the  nonlinear  prograimning  aspect 
of  the  maximum  principle.   Beyond  this,  it  might  also  prove  fruitful  to 
combine  the  latter  with  a  comparison  to  dynamic  programming. 

Finally,  the  discussion  of  numerical  instabilities  of  Chapter  VII 
indicates  the  possibility  that  the  respective  curvatures  of  the  initial 
guess  and  the  solution  determine  convergence.   Any  development  that 
could  indicate  convergence,  or  the  lack  thereof,  for  an  initial  guess 
required  by  an  indirect  solution  technique,  would  be  a  significant 
development.   One  possibility  is  a  study  of  the  sensitivity  of  the 
solution  increment  function,  £(t),  to  error  sources  such  as  initial 
guess,  step  size,  etc.,  for  a  case  with  a  known  solution.   Little  is 
known  about  this  subject. 


APPENDIXES 


APPENDIX  A 


HISTORICAL  DEVELOPMENTS 


The  dissertation  subject  matter  originated  from  two  papers: 
a  short  survey  of  structural  optimization  problems  by  Prager  and  Taylor 
(1968)  and  a  technical  note  by  Boykin  and  Sierakowski  (1972) . 
A  survey  of  the  literature  began  with  the  former  and  revealed  that  what 
were  presumed  to  be  recent  developments  were  in  fact  two  hundred  years 
old.   Besides  engendering  a  real  sense  of  humility,  my  discovery  aroused 
a  fear  that  some  fundamental  part  of  my  work  might  be  a  duplication  of 
a  much  earlier  study.   I  also  found  the  historical  development  to  be 
quite  intriguing;  for  that  reason  this  appendix  has  been  included  to 
give  the  reader  a  brief  history  of  the  separate  development  of  mechanics 
and  the  calculus  of  variations.   Neither  is  intended  to  be  comprehensive; 
however,  the  reader  may  enjoy  being  able  to  associate  some  recognizable 
equation  or  method  with  a  specific  person  and  period  of  time. 

The  appendix  consists  of  three  sections,  where  the  first  is  an 
anecdotal  discussion  of  the  major  early  contributors  to  the  calculus 
of  variations  and  a  table  listing  major  developments  in  chronological 
order.   This  is  follox^7ed  by  a  similar  account  for  mechanics,  also  having 
a  table  of  chronologically  ordered  developments.   Concluding  the  appen- 
dix is  a  short  biographical  table  that  lists  the  life  span  and  national- 
ity of  most  of  the  people  mentioned. 
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Newton  posed  a  problem  in  Book  II  of  the  Principia 
vJhich  requires  the  techniques  of  the  calculus  of  variations  for  solution. 
His  goal  was  to  find  the  solid  of  revolution  which  has  the  minimum 
resistance  in  axial  flow.   No  apparent  significance  was  attached  to  the 
problem  of  either  Newton  or  his  contemporaries. 

The  calculus  of  variations  actually  begins  with  two  Swiss 
mathematicians,  James  Bernoulli  and  his  younger  brother  John.   Until 
1690  the  latter  was  a  student  of  James,  but  soon  became  a  rival. 
In  June,  1696,  John  Bernoulli  posed  the  brachistochrone  problem  in 
Acta  Eruditorum: 


A  New  Problem,  to  the  Solution  of  I'/hich 
Mathematicians  Are  Invited 

Given  two  points  A  and  B  in  a  vertical  plane,  to  find  for 
the  movable  (particle)  M,  the  path  AI'B  descending  along  which 
by  its  own  gravity,  and  beginning  to  be  urged  from  the  point  A, 
it  may  in  the  shortest  time  reach  the  other  point  B. 


According  to  various  sources,  solutions  were  offered  by  Newton, 
Leibnetz,  de  I'Hopital,  and  the  Bernoullis.   At  Leibnetz'  request 
the  solution  was  withheld  to  encourage  others  to  consider  the  problem. 
In  January,  1697,  the  problem  was  reannounced.   In  the  following  May 
issue  of  Acta  Eruditorum,  the  solutions  of  the  Bernoullis  and  the 
Karquis  de  I'Hopital  were  published.   John's  is  alleged  to  be  the  most 
readable  but  since  it  is  solved  as  an  analog  to  an  optics  problem,  the 
method  is  restricted  to  a  small  class  of  problems.   The  solution  of 
James  is  quite  geometrical,  treating  the  curve  as  a  polygon  with  sides 
of  infinitesimal  length.   It  was  also  assumed  that  whatever  optimal 
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property  the  entire  curve  possessed,  was  also  possessed  by  each  part. 

The  solution  of  de  I'Hopital  was  given,  but  without  proof. 

At  the  conclusion  of  his  paper,  James  posed  two  more  difficult 

problems: 

First:   find  the  curve  of  quickest  descent  from  a  fixed  point  A 
to  a  given  vertical  line.   This  version  of  the  brachistochrone 
problem  involves  conditions  for  minimizing  a  definite  integral 
where  variations  at  the  limits  are  permitted.   Such  conditions 
lead  to  natural  boundary  conditions,  or  transversality  condi- 
tions, and  were  not  obtained  until  some  years  later  by  Lagrange. 

Second:   among  all  curves  with  given  length  and  base,  find  the 
the  curve  such  that  a  second  curve,  whose  coordinates  are 
some  function  of  the  first,  contains  a  maximum  or  minimum  area. 
This  subsidiary  condition  is  called  an  isoperimetrical  constraint. 

Besides  inviting  all  mathematicians  to  attempt  solutions,  James  specif- 
ically offered  his  brother  John  a  prize  of  fifty  ducats  for  the  correct 
solution — according  to  Bliss  (1925,  p.  12).   John  claimed  the  award  but 
James  refused,  showing  that  the  assumption  of  a  uniform  optimality  prop- 
erty for  each  infinitesimal  element  was  not  valid.   James  published  the 
correct  solutions  in  the  May,  1701,  issue  of  Acta  Eruditorum;  to  include 
the  effect  of  the  constraint,  he  used  three  adjacent  polygon  elements  of 
Infinitesimal  length  instead  of  two  that  sufficed  for  the  brachistochrone, 
and  no  fewer  than  four  solution  techniques: 

(1)  Equilibrium  of  forces  on  an  element 

(2)  Equilibrium  of  moments  on  an  element 

(3)  Principle  of  virtual  work 

(4)  Principle  of  minimum  potential  energy 

I'Jhether  or  not  the  prize  existed  in  immaterial;  it  is  true  the  problem 
embroiled  the  brothers  in  a  bitter  quarrel,  and  together  with  the  brach- 
istochrone problem  marks  the  beginning  of  the  calculus  of  variations. 
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All  of  this  Xv7ork  by  the  Bernoullis  was  based  upon  complicated 
geometrical  results.   A  student  of  theirs,  Leonhard  Euler,  expanded  and 
generalized  their  initial  developments.   His  first  major  contribution 
to  the  calculus  of  variations  was  the  differential  equation  bearing  his 
name.   Euler  suggested  that  for  continuous  functions  a  definite  integral 
vanishes  identically  if  the  integrand  vanishes.   VJhen  the  integrand  is 
the  product  of  some  function  M(x,x,t)  and  an  arbitrary  variation   6x(t) 
the  function  must  vanish  if  the  variation  is  to  be  arbitrary.   Euler's 
equation  is   M(x,x,t)  =  0;  up  to  the  first  part  of  this  century  this 
vanishing  integrand  was  referred  to  as  "Euler's  Fundamental  Lemma" — 
"Euler's  equation"  is  the  more  recent  appellation.   He  also  was  the  first 
to  indicate  the  two  general  classes  of  optimal  solutions  obtained  from 
his  technique:  absolute  and  relative.   After  reading  a  memoir  by  Lagrange 
which  introduced  the  symbol  "5,"  Euler  adopted  its  use  for  his  o\-m. 
studies  and  named  it  in  an  article  entitled  "Elementia  Calculi  Varia- 
tionum."   Other  investigations  related  to  this  include  combining  the  cal- 
culus techniques  of  Newton  and  Leibnetz  via  differential  elements,  and 
what  is  apparently  the  earliest  study  of  integrability.   Euler's  list  of 
886  books  and  articles  also  includes  the  introduction  of  "E"  to  denote 
summation,  and  an  intuitive  description  of  continuous  functions — all  of 
which  are  necessary  elements  of  the  calculus  of  variations. 

As  mentioned  in  the  previous  paragraph,  Lagrange  introduced  the 
concept  and  symbol  "6";  using  it  enabled  him  to  derive  Euler's  equation 
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without  resorting  to  the  infinitely-sided  polygon  and  the  accompanying 
complicated  geometry.   It  also  facilitated  the  evaluation  of  problems 
having  definite  integrals  with  unspecified  conditions  at  the  limits — 
these  were  called  "the  definite  equations,"  and  later  "terms  at  the 
limits,"  but  are  currently  referred  to  as  "natural  boundary  conditions." 
This  represented  a  generalization  on  the  type  of  problem  originally 
posed  by  James  Bernoulli.   Lagrange  introduced  yet  another  generaliza- 
tion by  using  "indeterminate  multipliers"  to  enforce  implicit  constraints, 
e.g.,  differential  equations.   The  final  major  contribution  by  Lagrange 
to  be  cited  here  was  the  extension  to  problems  having  double  integrals 
with  fixed  conditions  at  the  limits. 

Following  a  study  of  the  major  works  of  Newton,  Laplace,  and 
Lagrange,  William  Rowan  Hamilton  applied  the  analytical  method  of 
Lagrange  to  his  o\.m  study  of  optics.   In  the  "Theory  of  Systems  of  Rays" 
(1828)  and  its  supplements,  Hamilton  deduced  all  optical  properties 
from  a  single  "Characteristic  Function"  using  the  "Law  of  Varying  Action," 
a  generalization  of  the  "Law  of  Least  Action"  due  to  P .  L.  M.  de  Maupertis. 
When  applied  to  dynamical  systems  in  "On  General  Methods  of  Dynamics" 
(1834)  this  method  is  recognized  as  Hamilton's  Principle.   The  canonical 
equations  of  motion  derived  from  this  method  are  called  Hamilton's  equa- 
tions, and  were  used  by  Jacobi  to  later  derive  the  Hamilton-Jacobi  equa- 
tion. 

Euler  derived  the  first  necessary  condition  for  an  optimal  solu- 
tion in  the  calculus  of  variations:   in  the  manner  developed  by  Lagrange, 
the  first  variation  must  vanish.   But  this  may  be  satisfied  by  either 
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a  minimum,  maximum,  or  an  "inflection  point."  Legendre  made  the  first 
attempt  to  distinguish  what  the  case  may  be  by  investigating  the  second 
variation  in  1786.   Using  a  transformation  of  variables,  he  arrived  at 
what  was  believed  to  be  a  sufficient  condition  (for  either  a  maximum  or 
a  minimum)  for  fifty  years.   However,  in  1837  Jacobi  discovered  that 
Legendre's  condition  failed  in  some  cases.   The  "conjugate  points"  asso- 
ciated with  these  cases  were  found  from  derivatives  of  the  solution  to 
Euler's  equation  with  respect  to  the  constants  of  integration.  As  a 
result,  Legendre's  condition  was  determined  to  be  only  necessary,  along 
with  the  newer  Jacobi' s  condition,  and  not  sufficient  as  originally 
thought. 

Derivation  of  a  true  sufficient  condition  was  done  by  Weiers trass 
in  1879.   It  occurred  as  a  result  of  his  re-examination  of  earlier  works, 
formulating  the  problems  very  precisely.   Both  his  condition  and  his 
"corner  condition"  for  a  discontinuous  solution  are  stated  in  terms  of 
the  "Weierstrass  function." 

Further  elaborations,  involving  geodesies  and  the  existence  of 
integrals,  are  attributed  to  Kneser  and  Hilbert  about  1900.   Twenty  years 
later  Bliss  applied  the  "adjoint  system"  to  the  solution  of  a  variational 
problem.   The  "maximum  principle"  was  suggested  in  1937  by  Valentine  but 
is  generally  attributed  to  Pontryagin  in  1956. 

All  of  these  basic  developments  are  presented  in  Table  A.l, 
together  with  less  pertinent  details.   According  to  Bliss  (1925,  p.  181) 
the  classical  memoirs  of  the  Bernoullis,  Euler,  Lagrange,  Legendre,  and 
Jacobi  are  contained  in  Ostwald's  "Klassiker  der  exacten  Wissenschaf ten," 
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TABLE  A.l 
MAJOR  DEVELOPMENTS  IN  THE  CALCULUS  OF  VARIATIONS 


1686   Newton 


first  to  pose  a  problem  of  a  type  requir- 
ing the  calculus  of  variations;  not 
treated  as  significant  at  the  time. 


1696    John  Bernoulli 


1701   James  Bernoulli 


1722   Euler 


1740   Euler 


1744   Euler 


1756   Euler 


*   Lagrange 

■'•   Lagrange 


posed  the  brachistochrone  problem.   Solu- 
tions were  obtained  by  James  Bernoulli, 
Leibnetz,  and  Newton,  as  well  as  by  John 
Bernoulli  and  de  I'Hopital.   Accompany- 
ing his  solution  were  two  additional 
problems  proposed  by  James  Bernoulli, 
one  involving  natural  boundary  conditions, 
another  with  an  isoperimetrical  constraint. 

solved  the  isoperimetrical  constraint 
problem  posed  in  his  solution  to  the 
brachistochrone . 

extended  the  work  of  the  Bernoulli  brothers 
to  several  classes  of  problems  with  var- 
ious constraints. 

introduced  the  technique  that  setting 
coefficient  of   6Y  in  integrand  to  zero 
satisfies  necessary  condition.   The 
result  is  called  Euler' s  Equation. 

divided  optimization  problems  into  two 
general  classes:  absolute  and  relative. 
This  is  the  first  systematic  exposition; 
all  preceding  work  represented  ad  hoc 
solutions  to  specific  problems. 

published  one  tract  giving  the  name 
"calculus  of  variations,"  and  another 
containing  the  first  study  of  conditions 
of  integrability. 

introduces  the  symbol   "6"  to  distinguish 
between  variation  and  derivative. 

deduced  Euler 's  Equation  without  using  adja- 
cent elements  of  an  infinitely-sided  polygon. 


"no  specific  dates  are  available;  two  memoirs  on  the  general  subject  from 
which  these  might  originate  were  published  in  1762  and  1770. 
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A 


* 


Lagrange 


Lagrange 


Lagrange 


1786   Legendre 


1810   Brunaci 


determined  the  method  to  treat  unspec- 
ified end  conditions;  the  results  are 
now  called  natural  boundary  conditions. 

introduced  enforcement  of  implicit  con- 
straints by  the  use  of  indeterminate 
multipliers. 

extended  the  calculus  of  variations  to 
include  double  integrals  with  fixed 
boundary  conditions. 

investigated  the  second  variation  of 
the  cost  function  to  find  a  criterion 
with  which  maxima  and  minima  could  be 
distinguished.   First  results  were  not 
really  a  sufficient  condition,  because  of 
non-unique  solutions  demonstrated  by 
Lagrange 

extended  Legendre 's  work  to  double 
integrals,  retaining  the  same  general 
flaw. 


1834   Hamilton 


developed  the  canonical  equations  and 
Hamilton's  Principle  based  upon  a  general- 
ization of  the  Law  of  Least  Action. 


1837   Jacobi 


1842    Sarrus 


studied  the  transformation  of  Legendre, 
discovering  how  to  determine  when  the 
method  failed.   The  resulting  conjugate 
points  are  determined  from  derivatives 
of  solutions  to  Euler's  equation  with 
respect  to  the  constants  of  integration. 

won  the  French  Academy  of  Science  prize 
for  mathematics  by  obtaining  the  natural 
boundary  conditions  for  a  general  multiple 
integral  having  variable  limits.   His 
technique  used  a  new  substitution  sign  to 
designate  a  particular  value  for  general 
variable. 


1844   Cauchy 


simplified  the  work  of  Sarrus  to  the  form 
used  currently. 


No  specific  dates  are  available;  two  memoirs  on  the  general  subject 
from  which  these  might  originate  were  published  in  1762  and  1770. 
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TABLE  A.l  (Concluded) 


1852 


1871 


Brioschi 


1858   Clebsch 


Todhunter 


was  the  first  to  apply  determinants  to 
the  investigation  of  second  order  terms. 

generalized  Jacobi's  Theorem  to  several 
dependent  variables,  Vv^ith  or  without 
connecting  equations,  and  for  multiple 
integrals. 

suggested  that  variations  might  be  of 
restricted  sign,  allowing  the  possibility 
of  discontinuous  solutions. 


1877   Erdmann 


derived  the  conditions  that  must  be 
satisfied  by  a  discontinuous  solution, 
called  the  Weierstrass-Erdmann  conditions. 


1879   VJeiers  trass 


helped  to  make  the  theory  of  the  calculus 
of  variations  more  precise,  in  general. 
Specifically,  he  found  a  new  necessary 
condition  in  terms  of  his  function 
E(X,Y,P,Y')  and  very  clearly  determined 
those  conditions  which  are  sufficient. 


1919   Bliss 


1937   Valentine 


was  perhaps  the  first  to  use  the  adjoint 
system  of  variational  equations. 

derived  the  earliest  equivalent  to  a 
strong  minimum  principle,  using  the 
Weierstrass  condition  and  slack  variables. 


1956   Pontryagin 


together  with  colleagues,  derived  the 
Pontryagin  maximum  principle. 
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Nos.  46  and  47.   Historical  developments  are  treated  by  Carll  (1881), 
Todhunter  and  Pearson  (1893),  Bolza  (1904),  and  Bliss  (1925,  1946). 
Biographical  data  are  found  in  Truesdell  (1968),  Bergamini  (1963), 
Ireland  (1962),  and  Asimov  (1960). 

Analytical  mechanics  became  a  viable  method  with  the  introduc- 
tion of  calculus  by  Leibnetz  and  Newton.   For  example,  Galileo  Galilei 
observed  in  1585  that  the  free  fall  of  various  objects  of  different  size 
and  weight  is  given  by  the  formula  Y  =  16t^.   Newton  is  reported  to  have 
differentiated  twice  to  discover  that  all  objects  fall  with  an  acceler- 
ation of  thirty-two  feet  per  second  squared.   Yet  there  were  differences 
betwefen  the  calculus  of  Newton  and  Leibnetz;  the  former  designated  a 
derivative  by  a  dot  above  the  variable,  and  the  latter  used  the  form 
"dY/dt"  to  indicate  the  derivative  of  "Y"  with  respect  to  "t."   Newton 
also  deduced  that  maxima/minima  occur  at  points  where  the  rate  of  change 
is  zero,  whereas  Leibnetz  insisted  such  points  occurred  where  the  tangent 
to  the  curve  has  zero  slope. 

To  show  that  the  work  of  these  two  was  equivalent  was  one  of 
the  many  accomplishments  of  Euler .   In  using  this  discovery  he 
demonstrated  that  much  of  what  was  then  considered  to  be  pure  physics 
was  transformed  quite  simply  into  problems  of  mathematics.   Such  work  led 
him  to  sufficiently  important  developments  that  specific  equations  bear 
his  name  in  the  separate  fields  of  the  calculus  of  variations,  differ- 
ential equations,  solid  mechanics,  and  rigid  body  dynamics.   According  to 
Truesdell  (1968),  Euler  is  the  major  source  of  rational  analysis.   Cer- 
tainly he  is  one  of  the  most  prolific,  having  published  886  books  and 
articles,  and  fathering  thirteen  children. 
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A  more  detailed  list  of  the  basic  developments  of  mechanics 
is  given  in  Table  A. 2;  a   few   biographical  data  are  presented  in 
Table  A. 3. 

TABLE  A. 2 
MAJOR  DEVELOPMENTS  IN  ANALYTICAL  MECHANICS 


1694   James  Bernoulli 


1704   James  Bernoulli 


-*   James  Bernoulli 


1704   James  Bernoulli 


Leibnetz 


1727   Euler 


published  the  very  first  paper  on  the 
mathematical  theory  of  elasticity. 

concluded  a  thirteen  year  study  which 
included  an  investigation  of  the 
catenary  (results  were  published  in  1744) . 
Problem  was  also  solved  by  John  Bernoulli 
using  a  theorem  concerning  the  center  of 
gravity  of  an  arc. 

declared  elasticity  is  generally  non- 
linear, with  linear  elasticity  a  special 
case.   He  also  derived  the  differential 
equation  for  the  elastica. 

first  characterized  a  material  with  stress 
as  a  function  of  strain  instead  of  the 
then  accepted  load  versus  displacement 
for  a  particular  specimen.  He  avoided 
the  E  of  linear  elasticity,  suggesting 
E   by  implication. 

performed  the  first  analysis  of  a  loaded 
beam,  assuming  that  some  fibres  were  in 
tension.   That  the  tension  varied  lin- 
early gave  the  result  that  the  moment 
was  proportional  to  the  second  area  mo- 
ment of  the  cross  section. 

derived  the  Bernoulli-Euler  equation  for 
bending  of  beams  while  a  student  of  John 
Bernoulli.   The  work  was  not  published 
until  1862. 


No  specific  dates  are  available  from  ordinary  sources;  items  are 
located  with  respect  to  related  entries  of  the  table. 
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TABLE  A. 2  (Continued) 


1729   Musschenbroek     . . .    performed  the  first  comprehensive  exper- 
iments on  the  strength  of  materials.   He 
observed  that  beams  and  columns  bend 
before  breaking  and  that  breaking  strength 
varies  inversely  with  the  square  of  the 
length  and  also  depends  upon  depth  and 
breadth. 

1736    Euler  ...    introduced  the  study  of  mechanics,  specif- 

ically the  concepts  point  mass,  acceler- 
ation, and  vectors  (which  he  called 
"geometrical  quantity"). 

1739  Daniel  Bernoulli   ...    developed  the  principle  of  minimum  poten- 

tial energy,  although  he  is  best  knoiTn 
for  his  work  in  fluid  mechanics. 

1740  _  Euler  ...    learned  how  to  solve  analytically,  linear 

ordinary  differential  equations  with  con- 
stant coefficients.   Prior  to  that  time 
only  series  solutions  were  available. 

*   Euler  ...    after  making  some  calculations  on  an 

actual  beam  deflecting  under  its  load, 
discovered  that  most  real  applications 
result  in  small  strain.   All  previous 
work  had  presumed  that  deformations  were 
finite. 

1743   Euler  ...    determined  that  the  equation  of  the 

elastic  curve  could  be  obtained  from  the 
minimum  potential  energy. 

1743   Euler  ...    derived  the  column  buckling  formula;  the 

results  were  compared  to  Musschenbroek 's 
experiments  in  1765. 

1743   John  Bernoulli     ...    obtained  the  first  differential  equation 

for  the  motion  of  a  system  by  studying 
the  dynamics  of  a  weighted  cord  using 
lumped  masses. 


No  specific  dates  are  available  from  ordinary  sources;  items  are 
located  with  respect  to  related  entries  of  the  table. 
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TABLE  A. 2   (Continued) 


1743   D'Alembert 


1744   Euler 


1747   Euler 


1750   Euler 


introduced  the  first  partial  differential 
equations.   He  also  was  the  first  to 
separate  constraint  forces  from  external 
forces,  but  never  stated  "D'Alembert 's 
Principle,"  which  was  given  later  by  Euler 
and  Lagrange. 

published  the  earliest  example  of  what 
came  to  be  called  the  "Newtonian  Method." 
Also  defined  were  linear  momentum  and 
kinetic  energy  derived  by  integration. 

was  the  first  to  derive  "Newton's  Equa- 
tions" for  a  discrete  system.   By  apply- 
ing them  to  each  part  of  the  system,  all 
the  required  equations  were  obtained. 
The  example  used  was  the  three-body 
problem  of  celestial  mechanics. 

extended  the  preceding  study  to  "mechan- 
ical systems  of  all  kinds,"  i.e.,  to  con- 
tinuous systems.   Examples  used  included 
rigid  body  motion  and  transverse  vibra- 
tions of  a  rod. 


1754   Riccati 


1760   Euler 


*   Lagrange 

1773   Lagrange 


suggested  that  the  work  done  to  deform 
an  elastic  body  was  recoverable  in  part. 

developed  an  entire  theory  of  rigid  body 
motion.   He  defined  center  of  mass  and 
distinguished  it  from  center  of  gravity. 
Additionally,  he  obtained  the  inertia 
tensor  and  described  its  elements  with 
respect  to  the  mass. 

determined  the  infinite  sequence  of  theo- 
retical buckling  loads. 

attempted  to  find  the  column  shape  which 
would  have  the  largest  buckling  load  for 
given  height  and  volume.   His  analysis  was 
incorrect;  this  results,  a  circular 
cylinder. 


No  specific  dates  are  available  from  ordinary  sources;  items  are 
located  with  respect  to  related  entries  of  the  table. 
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TABLE  A. 2   (Concluded) 


1778   Euler 


1788   Lagrange 


1877   Gauss 


determined  the  height  at  which  a  uniformly- 
heavy  column  buckles. 

published  one  of  the  most  important  works 
in  mechanics:   "Mechanique  Analitique." 
In  this  book  he  systematized  mechanics, 
using  the  calculus  of  variations  to 
derive  general  equations  which  apply  to  all 
problems  of  mechanics.   Generally  speak- 
ing, he  obviated  the  use  of  complicated 
geometry,  replacing  it  with  pure  algebra. 
More  specifically,  he  introduced  the 
misnomer  "D' Alembert's  Principle,"  gen- 
eralized coordinates,  and  Lagrange's 
equations. 

proved  the  fundamental  theorem  of  algebra 
in  his  doctoral  dissertation,  requiring 
the  introduction  of  complex  numbers.   He 
later  recognized  their  vector  character 
and  applied  them  to  the  solution  of 
complicated  mechanics  problems. 
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TABLE  A. 3 

MAJOR  CONTRIBUTORS  IN  THE  DEVELOPMENT  OF 
THE  CALCULUS  OF  VARIATIONS  AND  ANALYTICAL  MECHANICS 


NEWTON,  Sir  Isaac  (1642-1727) 

LEIBNETZ,  Baron  von  (1646-1716) 

Gottfried  Wilhelm 

BERNOULLI,  John  (Jean,        (1667-1748) 
Johan) 

BERNOULLI,  James  (Jacques,     (1654-1705) 
Jakob) 

MUSSCHENBROEK,  Pieter  van     (1692-1761) 


BERNOULLI,  Daniel  (1700-1782) 

EULER,  Leonhard  (1707-1783) 

ALEMBERT,  Jean  Le  Rond  D'  (1717-1783) 

LAGRANGE,  Comte  Joseph  (1736-1S13) 
Louis 

LEGENDRE,  Adrian  Marie  (1752-1833) 

GAUSS,  Karl  Friedrich  (1777-1855) 

HAMILTON,  Sir  William  (1805-1865) 
Rowan 

JACOBI,  Karl  Gustav  Jacob  (1804-1851) 

BRIOSCHI,  Francesco  (1824-1897) 

CLEBSCH,  Alfred  Rudolf  (1833-1872) 
Friedrich 

TODHUNTER,  Isaac  (1820-1884) 

WEIERSTPxASS,  Karl  Theodor  (1815-1897) 

BLISS,  Gilbert  Ames  (1876-1951) 

PONTRYAGIN,  Lev  S.  (1908-    ) 


English  mathematician 
German  mathematician 

Swiss  mathematician 

Swiss  mathematician 

Dutch  physicist  and 
mathematician 

Swiss  mathematician 

Swiss  mathematician 

French  astronomer  and 
mathematician 

French  mathematician 

French  mathematician 

German  astronomer  and 
mathematician 

Irish  m.athematician 

German  mathematician 
Italian  mathematician 
German  mathematician 

English  mathematician  and 
historian 

German  mathematician 

American  mathematician 

Russian  mathematician 


APPENDIX  B 
A  SIMPLE  PROOF  OF  THE  KUHN-TUCKER  THEOREM 

Consider  the  minimization  of  a  general  unspecified  functional 
with  respect  to  an  n-dimensional  vector  x>  subject  to  m  <  n  inequal- 
ity constraints. 


functional:  Min  [F(x)] 

inequality  constraints:   g . (x)  <  0  j  =  1, . . . ,m 


where 


—     12      n 


T 


Inequality  constraints  are  converted  to  equality  constraints  by  intro- 
ducing "slack"  variables   s  .  (^c)  defined  such  that 

gj(x)  +  s2(x)=  0     j  =  l,...,m         (B.l) 

or 


5^ 

'2 


(B.2) 


s^(x)  =  [-gj(x)]^   >  0 

Since  for  all  j  the  left-hand  side  of  equation  (B.l)  sums  to  zero, 

for  arbitrary  values  of  A. 

J 

m 

E   A.[g.(x)  +  s2(x)]  =  0  (B.3) 

j=l  J   J       J 

it  follows  that: 


m 


$(x,s,A)  =  F(x)  +  Z   A.[g.(x)  +  s2(x)]  =  F(x) 

j  =  l   -^   ^       J 
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Necessary  conditions  for  the  adjoined  cost  function  $  to  be 
a  minimum  are  derived  from  the  requirement  that  the  first  variation 
vanish. 

8$  =  OX H  OS   -7 h    OA   — T  =  0 

If  all  the  variations  are  to  be  independent,  then 

V  F(x)  +  A.  V   g.(x)  =  0        i  summed 

2A.  s . (x)  =0        i  =  1, . . . ,m 
J   J  - 

g. (x)  +  s^(x)  =  0        j  =  l,,..,m 

where  repeated  indices  in  a  product  imply  summation.   The  independence 
of  _s  and  _A  can  be  argued  heuristically.   Slack  variables  are  an  arti- 
ficial contrivance  used  to  force  a  desired  condition  without  altering 
the  original  problem.   Any  change  in  s . (x)  which  satisfies  the  given 
constraint  must  be  accompanied  by  a  change  in  g.(x)  such  that  the  net 
effect  on  0  is  zero;  thus  *(2<,s^,_A)  is  stationary  with  respect  to  s^. 
Furthermore,  since  the  sum  in  (B.3)  is  equal  to  zero  for  any  arbitrary 
values  of  A.,  any  arbitrary  variation  of  $  with  respect  to  A_  must  also 
be  zero. 

Given  the  definition  of  s . (x)  in  (B.2),  the  three  necessary 

J  - 

conditions  may  be  expressed  as 

V  [F(x)  +  A.g.(x)]  =  0  (B.4) 
X   -     J  J  - 

A.g.(x)   =  0  (B.5) 

g.(x)  <    0 
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The  first  condition  will  be  used  with  the  concept  of  a  feasible,  usable 
direction  to  derive  the  final  condition  remaining  before  the  Kuhn-Tucker 
condition  may  be  stated.   Any  direction  in  2£-space  is  feasible  if  an 
increment  Ax  in  that  direction  reduces  the  cost  function  F(x) .   Direc- 
tion d  is  said  to  also  be  usable  if  an  increment  A>:  in  that  direction 
causes  no  constraint  violation. 

Consider  first  the  concept  of  a  feasible  direction;  for  any 
direction  d  the  directional  derivative  is  defined  to  be 

dF 

^  =■-   d    '    V  F(x) 

ds       2S  ~" 

where 

Pll  =  1 

In  addition  to  this  restriction  on  the  norm  of  direction  3,  let  the 
increment  dx   be  prescribed  by 

dx  =  (ds)d  (B.6) 

and  the  vector  scalar  product  be  written  as  the  equivalent  inner  product 
such  that 

d  •  V  F  =  d""-  V  F  =  (V  F,  d) 
In  terms  of  an  inner  product  notation 

For  an  increment  in  direction  d  of  length  ds,  the  change  in  the  cost 
function  is 


dF  =  (V  F,  d)  ds 

X 
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By  the  linearity  property  of  the  inner  product 

dF  =  -  (-V  F,  d)  ds 

Hence  the  cost  T (x)    is  decreased  by  an  increment  in  any  direction  d 
which  satisfies 

(-V  F,  d)  >  0  (B.7) 

X 

Conversely,  a  direction  d  is  feasible  if  (-V  F,  d)  >  0,  i.e.,  the  angle 
between  d  and  the  direction  of  steepest  descent  is  acute. 

For  any  direction  to  be  usable,  an  increment  !\x   in  that  direc- 
tion must  cause  no  constraint  violation.   On  the  constraint  boundaries 
g.(x)  =  0  for  j  e  I  as  defined  in  Chapter  III.   V/ith  any  direction  again 
denoted  by  d,  the  change  in  g . (x)  with  respect  to  an  increment  dx  is 
defined  through  the  directional  derivative: 

dg. 


J-   = 


ds    ^  X  ''j 


(V   g.(x),d) 


where 


such  that 


lai!  =  1 


dg,  =  (V   g  ,  d)  ds 
J      ii   J 

Since  ds  denotes  a  vector  norm  by  way  of  the  definition  of  dx^  given  in 
equation  (B.6),  it  is  a  positive  quantity.   Hence  the  sense  of  change 
dg.  is  prescribed  by  the  inner  product  (V   g . ,  d)  for  which  there  are 
three  pjossibilities. 
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(i)   (V  g.,  d)  <  0  ->  dg.  <  0 

It  was  assumed  that  active  constraints  are  being  discussed,  for 
which  g .  (2£)  =  0.   The  constraints  require  all  g.(x)  to  be  non- 
positive.   An  increment  dx   in  the  direction  S  used  in  the  inner 
product  satisfies  the  constraint  with  the  inequality  sign. 
Because  the  direction  does  not  lead  to  a  constraint  violation 
it  is  by  definition  usable.   Actually,  for  such  a  direction  d, 
even  though  point  x  lies  on  the  boundary  g . (x)  =  0,  the  point 
(x  +  dx)  does  not. 

(ii)   (V  g.,  d)  =  0  ^  dg.  =  0 
2i  J  J 

Under  the  assumption  that  g.(x)  is  active  where  g . (x)  =  0,  when 
the  increment  d^  in  direction  d  causes  no  change  in  g.,  d  is 
tangent  to  the  constraint  surface.   With  x  +  "^J^  used  in  the 
inner  product,  the  point  (2c  +  dx^)  may  still  lie  on  the  constraint 
boundary.   With  such  an  increment  the  constraint  may  perhaps  not 
be  violated,  and  if  not,  the  direction  is  therefore  usable. 
Note  that  for  linear  constraints,  the  word  "may"  is  deleted  from 
the  previous  sentence.   It  is  included  because  nonlinear  con- 
straints require  an  argument  calling  for  infinitesimal  incre- 
ments; a  finite  increment  leads  to  a  point  that  is  not  on  the 
constraint  boundary.   A  finite  increment  from  a  point  on  a  non- 
linear constraint  boundary  results  in  that  constraint  either 
becoming  inactive  on  being  violated,  depending  upon  whether  the 
constraint  surface  is  concave  or  convex,  respectively. 
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(iii)   (V  g.,  d)  >  0  ^  dg.  >  0 

That  the  constraint  g.(x)    is  active  requires  g . (x)  =  0. 
If  an  increment  <ix  in  the  direction  d  used  in  the  inner 
product  results  in  a  positive  dg.,  the  constraint  has  been 
violated.   Therefore,  d  is  not  a  usable  direction. 


The  result  of  these  possibilities  is  that  for  any  point  x  on 
a  constraint  boundary,  for  a  direction  d  to  be  usable  it  must  satisfy 

(V  g.,  d)  <  0  j  e  I.  (B.8) 

If  the  set  of  all  directions  d  which  are  feasible  contains  no  usable 
directions  at  the  point  x   in  question,  then  x   =  ^^p^^   and  x   is  the 
constrained  minimum  point. 

Conditions  for  a  usable,  feasible  direction  given  by  (B.8)  and 
(B.7)  are  used  with  the  gradient  condition  (B.4)  to  prove  that  the 
Lagrangian  multipliers  of  (B.5)  must  be  non-negative  at  the  point  where 
X  is  a  constrained  optimum.   Condition  (B.4)  states  that  at  the  con- 
strained minimum  point,  the  negative  gradient  of  the  cost  function  lies 
in  a  subspace  defined  by  the  gradients  of  the  active  constraints,  i.e., 

-  V  F(x)  =  A.  V   g.(x) 
X         J   X   J  ~ 

In  the  following,  implied  summation  is  replaced  by  the  symbol  E,  and 

the  dependence  of  F  and  g.  upon  x   is  implied.   Resolving  -V  F  into  com- 

1  2£ 

ponents  along  the  V  g.  directions  via  the  inner  product  gives 
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m   (-V^F,  V^g  ) 


Thus 


-V  F  =  I     ^ =^  V  g. 


A.  =  = =-^   ,  J  e  I 

^      I|V  g.  II 


A 


This  condition  relating  the  Lagrangian  multipliers  to  an  inner 

product  is  used  x^7ith  the  conditions  necessary  for  the  existence  of  a 

constrained  minimum  to  show  that  A.  >  0.   At  a  constrained  optimal  point, 

J 

for  every  feasible  direction  d,  the  cost  decreases  and 
(-V  F,  d)  >  0  ^   d  is  feasible 

X 

the  constraints  are  violated  (otherxrise  not  optimal). 


(V  g.,  d)  >  0  <-     d  is  not  usable 

Assume  that  some  j,  A.  is  negative  and  a  feasible  but  not  usable  direc- 
tion S.   exists.   If  d  is  taken  to  be  coincident  with  -V  F,  this  direc- 

x^ 

tion  is  feasible  and  usable,  contrary  to  the  original  assumption. 
That  A.  >  0  is  proved  by  contradiction.   Mathematically  stated, 
Given:      (-V  F,  d)  >  0  for  all  d 

(V  g.,  d)  >  0  for  all  j  e  I 

X  J  A 

Assume:      A.  <  0  for  some  j  e  I. 

(i)      A.  <  0  ->  (-V  F,  V  g.)  <  0 

J  X      X  J 


213 


(ii)     Consider    d  =  --V  F 

X 


then 


(-V  F,  d)  =  (-V  F,  -V  F)  =   ll-V  F  l|2 

X  XX        "    X   " 


>  0  ^  d   is  feasible 


(iii)    But 


(V^g.,  d)  =  (V^g.,  -V^F) 


<  0  ^  d   is  usable 


by  virtue  of  (i) .   This  says  the  direction  of  steepest  descent 
is  a  feasible  direction  which  is  also  usable,  contrary  to  the 
given  condition. 


(iv)     Therefore, 


A.  >  0,      j  e  I,  Q.E.D. 

.1  A 


This  completes  the  proof  of  the  Kuhn-Tucker  conditions  for  the 
existence  of  constrained  minimum  of  Fix)    subject  to  g .  (>c)  <  0 
j  =  l,...,m.   They  are 


(i)  gjCxop^)  ^  0 


(ii)  A.  >  0 

J 


^j8j(x,p,)  =  0 


m 


(iii)        \H^^^)   +  _Z   A   V^g  (x^p^)  =  0 


APPENDIX  C 
COMPUTER  SUBROUTINE  LISTINGS 

In  order  to  use  the  IBM  SHARE  quasilinearization  program,  three 
user-supplied  subroutines  must  be  combined  vvfith  three  SHARE  subroutines. 
The  latter  three  are  written  in  terms  of  a  general  problem  described  in 
section  3.5  and  do  not  change  from  one  problem  to  another.   Their  names 
and  functions  are: 

QUASI   -  a  quasilinearization  algorithm 

RECIP   -  a  m.atrix  inversion  subroutine 

OUTPUT  -  a  subroutine  that  prints  out  different,  specified 
combinations  of  problem  variables 

In  addition  to  these  IBM  subroutines,  the  user  must  provide  three  sub- 
routines which  convert  a  specific  problem  into  the  general  formulation 
required  by  the  above  subroutines.   Their  names  and  functions  are: 

MMN   -  reads  data,  establishes  initial  guess,  calls  QUASI, 
and  terminates  execution  xd.th  additional  output  if 
desired 

DIFEQ  -  J.s  called  by  QUASI  to  evaluate  the  differential 

equation  and  linearization  gradients  at  every  point 
for  each  iteration 

CORRC  -  corrects  specified  boundary  conditions  after  each 
iteration 

These  six  subroutines  are  combined  into  a  single  computer  program  for 

solving  the  given  TPBVP. 


214 


215 


Listings  of  nine  subroutines  are  contained  in  this  appendix. 
The  first  six  represent  a  program  for  solving  the  minimum  deflection 
beam  problem  subject  to  stress  constraints.   The  order  of  presenta- 
tion is: 

MAIN 

DIFEQ 

CORRC 

QUASI 

RECIP 

OUTPUT 

Following  these  are  the  three  use-supplied  subroutines  required  to 
solve  the  maximum  buckling  load  problem  subject  to  only  geometric  and 
material  bounds.   They  are  presented  in  the  following  order: 

MAIN 

DIFEQ 

COPJIC 
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